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ABSTRACT 

We present a theoretical framework for formal study of systematic effects in Supernovae Type la (SN la) 
that utilizes two-dimensional simulations to implement a form of the deflagration-detonation transition (DDT) 
explosion scenario. The framework is developed from a randomized initial condition that leads to a sample of 
simulated SN la whose ^^Ni masses have a similar average and range to those observed, and have many other 
modestly realistic features such the velocity extent of intermediate mass elements. The intended purpose is 
to enable statistically well-defined studies of both physical and theoretical parameters of the SN la explosion 
simulation. We present here a thorough description of the outcome of the SN la explosions produced by our 
current simulations. 

A first application of this framework is utilized to study the dependence of the SN la on the ^^Ne content, 
which is known to be directly influenced by the progenitor stellar population's metallicity. Our study is very 
specifically tailored to measure how the ^^Ne content influences the competition between the rise of plumes of 
burned material and the expansion of the star before these plumes reach DDT conditions. This influence arises 
from the dependence of the energy release, progenitor structure and laminar flame speed on ^^Ne content. For 
this study, we explore these three effects for a fixed carbon content and DDT density. By setting the density at 
which nucleosynthesis takes place during the detonation phase of the explosion, the competition between plume 
rise and stellar expansion controls the amount of material in nuclear statistical equilibrium (NSE) and therefore 
^^Ni produced. Of particular interest is how this influence of ^^Ne content compares to the direct modification 
of the ^''Ni mass via the inherent neutron excess as discussed by Timmes, Brown & Truran (2003). Although 
the outcome following from any particular ignition condition can change dramatically with ^^Ne content, with 
a sample of 20 ignition conditions we find that the systematic change in the expansion of the star prior to 
detonation is not large enough to compete with the dependence discussed by Timmes, Brown & Truran (2003). 
In fact, our results show no statistically significant dependence of the pre-detonation expansion on ^^Ne content, 
pointing to the morphology of the ignition condition as being the dominant dynamical driver of the ^^Ni yield 
of the explosion. However, variations in the DDT density, which were specifically excluded here, are also 
expected to be important and to depend systematically on ~^Ne content. 

Subject headings: hydrodynamics — nuclear reactions, nucleosynthesis, abundances — supernovae: general 
— white dwarfs 



L EVITRODUCTION 

Type la supernovae (SN la) are bright stellar explosions 
spectroscopically distinguished by strong silicon features near 
maximum light and a lack of hydro gen features ( Filippenko] 
T997| [Hiilebrandt & Niemeyer||200 0). Motivated bro adly by 
the importance of SN la for cosmological studies ( Phillips | 
[1993] [Riess et al.'T996'; Alb recht et a l. 2006), contempo- 
rary observational campaigns are gathering information about 
SN la a t an unprecedented ra t e (BAOSS,|Li et al.| 1996 2001 ; 
LOSS, 'Treffer s et al] [19971 |L i et al.' ioOl'; SNLS, .Astier| 
erar|r2006i CSP, [Hamuy et al. 2006; Nearby SN Factory, 
Copin et al. 2006"; Skymapper, Keller et al. 2 007] ESSENCE 
Wood-Vasey et al.,,2007 ; STRESS, Botticella et al.| [SOPS' 



SDSS-II, [Hdtzman et al.|2008[ SXDS, [Furusawa et al.|2008l 
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'Totani et al. 
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|2009| Pan-STARRS; the La Silla SN Search) and the future 
pro mises even more (the Dark Energy Survey, LSST, JDEM, 
see Howell et al.||2009a I. Interesting systematics have been 
discovered and continue to be better characterized. All SN la 
appear to burn a similar amount of total material, but can 
differ widely in the amount of ^''Ni produced, an effect that 
is closely correlated with their brightness and decline time 
( [Mazzali et al.|2007] l. Observations indicate that there are two 
populations of SN la that differ in the elapsed time between 
star formation in the host galaxy and the explosion ( Scanna- 
[pieco & Bildsten 2005; Mannucci et al. 20061. These two 



populations also appear to have a different average brightness, 
with SN la occurring in active star-forming galaxies appear- 
ing brighter ( Howell et al.|[200"7 i. Recent observations find 
correlations suggesting that SNe la in galaxies whose popula- 
tions have a characteristic age greater than 5 Gyr are 1 mag 
fainter at maximum luminosity than those found in galaxies 
with younger populations, while progenitor metal abundance 
has a weaker influence on peak luminosity (I Gallagher et al.| 
[2008| l. The color properties of observed SNe la have be- 
come an important prospective source of systematic error for 
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cosmic measurements ( [Hicken et al. 112009b!). Alongside the 
finding that the scatt er from th e Hubble law depends on host 
galaxy metallicity ( jGallagher et al.|2008j) , this creates concern 
that corrections for extinction have become entangled with 
uncharacterized intrinsic color variation. It is in this context 
of extensive empirical searches for systematics that our broad 
but still incomplete theoretical understanding of these events 
must be pushed as far as possible, to attempt to predict and 
understand systematic trends. 

By invoking a straightforward counting argument based on 
the fact that the number of protons and neutrons are approx- 
imately conserved in the explosion, [Timmes et al.| ( |2003[ ) ar- 
gued that there should be a fairly robust metallicity effect on 
the average ^^Ni yield of SNe la, and therefore, potentially, 
their brightness. Motivated by the simplicity of this effect and 
the important implications for cosmological usage of SNe la, 
significant effort has gone into measuring such a metallicity 
dependence observationally, but a clear effect has proven elu- 
sive ([Gallagher et al. l2005 ; Gallagh ereFal.|2008| [Howell et al.| 
pOOQbi Constraining metallicity dependence alone is chal- 
lenging for two reasons beyond the fact that the effect does 
not appear to be large. It is fairly difficult to measure accu- 
rate metallicities for the parent stellar population, and even so 
there are strong systematic problems with such measurements 
due to the mass-metallicity relationship within galaxies (Gal- 
|lagher et"al. 200 8 ). Additionally there is an apparently much 
stronger dependence of mean brightness of SN la on the age 
of the parent stellar population ( Gallagher et al. 2008 Howell 
let al.|2009b^ . 

The most recent observational work leaves the situation still 
murky. Gallagher et al. ( 2008| were unsuccessful at finding 
a metallicity dependence of the average brightness, but did 
find a such a dependence of the Hubble residual (scatter from 
the Hubble law) obtained from light-curve fitting. In contrast, 
Howell et al.| (j2009b ) found a slight dependence of average 



brightness on metallicity and no similar Hubble residual is 
sue using a different light-curve fitting method. However, the 
effect is difficult to separate from the dependence on mean 
stellar age. Given such uncertainty, it is important to con- 
tinue to evaluate the influence other aspects of the explosion 
might have on the effect of the simple overall neutron excess 
outlined by [Timmes et al. ( [2003| l. The importance is further 
stressed because it is clear that other effects, including the in- 
trinsic random variation of otherwise similar SNe la, can be 
even larger in magnitude. 

The systematic effects of metallicity on the SN la outcome 
have been the topic of a number of theoretical studies with 
a variety of explosion models. Several of these were accom- 
pUshed in one dimension ([Hoflich et al.|1998 Iwamoto et al. 



T999||Hoflich et al.|2000|[Domfnguez et al.|2001| l, and there 
fore bear revisiting with multi-dimensional models. In addi- 
tion to changing the overall yield of ^^Ni, the initial neutron 
excess, as set by the metallicity, is important for the compo- 
sition of intermediate layers that act as an important opac- 
ity source during the photospheric phase of the supernova 
( Dommguez et al.[2001[ ). The challenges of understanding the 
SN la phenomenon with multi-dimensional numerical mod- 
els have since somewhat overshadowed this type of study of 
population systematics. For example, one multi-dimensional 
study by Ropke et al. (20061 studied systematics using a less 
parameterized supernova model, but as a result was forced 
to treat marginally successful deflagration models that do not 
produce realistic explosions. Additionally, the neutron excess 
was treated only in post-processing, thus excluding sensitivity 



to dynamical effects like those studied here. 

The number of possible systematic parameters to consider, 
along with the intrinsic scatter expected from implementation 
of realistic turbulence characteristics, leads us to develop a 
theoretical framework for formally evaluating systematic ef- 
fects and their possible statistical significance. Systematic ef- 
fects that require study include both physical ones such com- 
position and ignition density and purely theoretical ones such 
as detonation condition and flame model parameters. In addi- 
tion to the deflagration-detonation transition (DDT) explosion 
model itself, the basic component of this framework is an ini- 
tial condition that defines a theoretical population or ensemble 
from which we can draw a sample of supernovae and study 
how the sample as a whole responds to parameter changes, 
giving a more complete and statistically quantifiable picture 
of their systematic impact. 

The ultimate goal of the study of SN la systematics is to 
understand the dependence of properties of a SN la popula- 
tion on characteristics of the parent stellar population. Among 
other benefits, this understanding would allow known charac- 
teristics of the stellar populations of galaxies, such as their 
cosmic history, to be utilized in understanding systematic 
trends that may appear in SN la, or to better understand the 
chemical evolution of galaxies. Stellar populations are most 
basically characterized by their age and metallicity, or more 
realistically, some mixture of or distribution of these parame- 
ters. There are likely to be several other important secondary 
parameters such as binarity, or environment (e.g. ionizing 
background) that may have enough influence on a stellar pop- 
ulation to be reflected in its SN la population. The domi- 
nance of stellar population age and metallicity (|Scannapieco| 
& Bildsten 2005 ; Ma nnucci et a l. 2006; Gaflagher et al. 2()08[ 
iHowell et al.„2009 b), or, alternatively, host galaxy morphol- 
ogy or color, in observational investigations to date indicates 
that such secondary factors are probably small by compari- 
son. However, until the age dependence is fully understood, 
it is important to be mindful of possible additional contribu- 
tions. 

Our aim in the present paper is far more modest. We be- 
gin by considering the dependence of the supernova on the 
composition of the exploding WD, and refine this exploration 
to only particular dynamical aspects in order to allow a tar- 
geted, conclusive result. The question that we seek to answer 
is: does the change in the dynamics of the explosion due to a 
different ^^Ne abundance in the progenitor introduce a signif- 
icant systematic effect in addition to the neutron excess dis- 
cussed by [Timmes et al.| (|2003)? This highly targeted scope 
is motivated by a desire give a clear and extensive description 
of our framework. This is the first time that we are applying 
two-dimensional DDT simulations as a method for generating 
a semi-realistic theoretical sample of supernovae and studying 
the systematics of that sample. 

The DDT model proved to be one of the most success- 
ful of the one-d imensional SN la models (e.g. [Hoflich &] 
Khokhlov|1996 1. However, it was never satisfying that both 



the deflagration velocity and the DDT transition density were 
treated essentially as free parameters. The hope of multi- 
dimensional models is that burning propagation during the 
deflagration phase, which was necessarily parameterized in 
one dimension, can be calculated directly. This would re- 
move another free parameter and lead to more reliable mod- 
els. Unfortunately the manifestation of buoyancy instabili- 
ties in multi-dimensional models became a serious challenge. 
Even modest asymmetries in the initial conditions of the de- 
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flagration led to far too little expansion of the star by the 
time that a traditional DDT would occur (|Niemeyer et al.| 
[19961 ICalder et al.|2003|[2004||Livne et al.|2005» . ITiis can 
be ameliorated with particular choices of ignition condition, 
allowing the main desirable features of the 1-d models to 
also be obtained from multi-dimensiona l models (Golombek 
& Niemeyer||2005[ [Gamezo et al ||2005[ [Ropke & Niemeyer 



20071 [Bravo & Garcia-Senz|2008| l. 



The degree to which such a symmetric deflagration is phys- 
ical is still a matter of some debate. The distribution of flame 
ignition regions in time and space remains fairly uncertain 
( [Woosley~et al. 2004; Rop ke et al.|20 06 ) as well as the degree 
to which turbulence-flame interaction after ignition can influ- 
ence the subsequent spread of the flame ( Ropke et al.[[2007 



[Jordan et al.|2008| . When placed alongside the physical un 
certainty as to whether or not a transition to detonation can 
actually occur (Nie meyer[[ 1 999| l, in a certain light the DDT 
scenario may simply be too contrived. However, so far the 
evidence is not sufficient to disprove it, and it continues to 
provide one of the best prospects for matching the observa- 
tions of the typical SN la. 

We begin in Section [2[ by discussing the variety of physi- 
cal effects through which composition can lead to systematic 
variation of SN la properties. This forms the context for how 
and why we limit the scope of this first study to dynamic ef- 
fects only. Following this, in Section [3] some improvements 
to the numerical model of the explosion are discussed and the 
extensions necessary to treat detonation-flame interaction and 
the presence of ^^Ne are presented. The burning model, flame 
speed, and mesh refinement are each discussed. In Section[4l 
we present the ignition condition on which our ensemble of 
SNe la is built and use it to construct a framework for eval- 
uating the significance of systematic effects. The varieties 
of yield arising in the theoretical SN la population are dis- 
cussed along with how the initial condition appears to control 
the outcome. A metric for measuring the expansion prior to 
DDT is introduced and calibrated to reflect the mass of nu- 
clear statistical equilibrium (NSE) material synthesized. Sec- 
tion[5]presents a detailed description of the outcome of a rep- 
resentative two-dimensional simulation and how it generally 
compares with observations. Some points of the implemen- 
tation of the explosion model that necessitate referring to the 
details of the explosion are also described here. Finally, in 
Section [6] the framework is applied to measure the effects of 
^-Ne content on the expansion of the star at the time of the 
DDT. We then summarize our conclusions and discuss future 
work in Section [71 

2. SYSTEMATIC INFLUENCES OF COMPOSITION 

We begin with an overview of the physical effects via which 
composition can influence the process and outcome of a DDT 
explosion. The most important material constituents of the 
WD are '^C, '^O, and ^^Ne, and we will frame our discussion 
in terms of these. Composition can influence the explosion 
through changes in several physical processes and properties 
involved in the explosion. These include ignition density, 
DDT density, energy release, flame speed, WD structure, and 
neutron excess. Detailed treatment of the first two of these, 
ignition density and DDT density, will be deferred to future 
work. The others will all be treated here, with the last, neu- 
tron excess, taken as the baseline effect to which others should 
be compared. Several of these effects involve significant un- 
certainties, and it is useful to highlight each in turn. 

The compositional structure of the carbon-oxygen WD that 



explodes is determined principally by the post-main sequence 
evolution of the star of which it is a remnant ( Domfnguez et al.| 
2001| l. The inner several tenths of a Mq are formed during the 



convective core helium burning phase, and the layers outside 
this in s hell burning on the asymptotic giant branch ("St ranieroj 
et al.|2003, and references therein). During helium burning. 



the C, N, and O w hich was preserit in th e initial star is trans- 
formed into ^^Ne ('Timmes et al. 2003 i, leading to a direct 
dependence on metallicity of the parent stellar population. In 
addition, ^^Ne is formed in the pre-explosion convective car- 
bon burning core at a comparable abundance (,Piro & Bildsten] 
|2008t|Chamulak et al.|2007l l. 

Ignition density - The ignition density characterizes when, 
as the result of accretion, the WD core begins runaway heat- 
ing due to carbon fusion outpacing neutrino losses (Nomotoj 
,1982 ; Woosley & Weaver|1986 ). The central density at flame 
ignition, when convection is insufficient to spread the heating 



throughout the core, is slightly less (see e.g. Piro & Bildsten 
[2008| l, and is also often called the ignition density, the mean- 
ing usually being discernible from context. The flame ignition 
density is generally near ~ 2 x 10^ g cm"^ for successful mod- 
els of SN la (e.g. [Nomoto et al":i[T984l [H oflich & Khokhlov] 
[1996) 1. As we do here, [Hoflich et al.l ( |1998, ,2000 ) adopted a 
single value for the ignition density for their studies of com- 
position dependence. Although we exclude it here, the varia- 
tion of ignition density is expected to be a significant effect. 
Along with the energetic variations due to the carbon con- 
tent discussed below, it is likely an important contribution to 
the observed dependence on stellar population age, which has 
proven to be stronger than any metallicity dependence ( |Gal-| 



lagher et al..2008 ; Howell et al.,200 9b). 

The precise value of the ignition density is sensitive to sev- 
eral factors, each of which has remarkable uncertainties. One 
of the reasons we will simply fix the ignition density for this 
study is the variety of uncertainties which must be consid- 
ered if it is varied. The energy generation rate depends on the 
'^C abundance in the core, as set by the evolution of the star 
that formed the WD. Authoritative calculation of this abun- 
dance remains elusive due to its dependence on the details 
of con vection during late core helium burning (e.g. [Straniero| 



et al. 2003). The core temperature of the WD is also im- 



portant, and so the thermal history of the core, notably the 
accretion rate and possibly properties of the helium flashes 
( [Nomoto et al.[[T984| . Either carbon composition or thermal 
state could lead to metallicity dependencies that are closely 
involved with both the evolution of the parent stellar popula- 
tion and the still very poorly understood ( [Branch et al. [19951 ) 
process of progenitor system formation. There are also un- 
certainties in the screening enhancement of nuclear reactions 
at these high densities ( [Gasques et al.[|2005 [Yakovlev et al.| 
[2006| ), and in the reaction rates themselves. In particular, the 
C+ C reaction cross-section is not experimentally deter- 
mined down to the low center- of-mass energies relevant for 
ignition in white dwarfs. There is evidence, from heavy-ion 
fusion reactions, of "hinderance" — a suppression of the astro- 
physical S-factor — a t sub-barrier energies Oiang et al.||2007 
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Gasques et al. 2007). In the case of '^C+'^C, however, reso 



nances are predicted in the energy range of interest (Michaud] 
& Vogt|1972[ Perez-Torres et al.|2006[ ), which could signifi- 



cantly incre ase the cro ss section by as much as two orders of 
magnitude ^Cooper et al. 2009 ) at temperatures ^ 5 x 10^ K. 

DDT density - We make the presumption that there is a 
unique characteristic density at which there is a transition 
from deflagration to detonation. In future work, this will be 
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extended to comparison of flame width and turbulent state 



(INiemeyer & WoosleyjlQgTllKhokhlov et al.|1997||Golombek| 



& Niemeyer||2005| l, giving a more physical transition point. 
While the details of this transition remain difficult to fully 
quantify ( Aspden et al. 2008 ; Woo sley et al.| 2008 ), it is very 
likely that it will have direct dependencies on composition 
because both the reactivity and the flame width depend on 
both the '^C and ^^Ne abundances. Note that in this case, 
the important composition is that in the outer portions of the 
star. In addition to having a higher '^C abundance than the 
center (e.g. Hoflich et al.n2000 ) because it is the product of 
shell burning, this region will have a lower ^^Ne abundance 
because it will not have participated in the core conve ction 



during which ^^Ne is enhanced in the convection zone (Piro 
|& Bildsten|2008) IChamulak et al.|2 008 ). Given the outcome 
of this study, showing that dynamical factors have little im- 
pact on the NSE yield, the influence on the DDT density is 
expected to be the principal avenue, beyond inherited neu- 
tron exces s, for dependence on metallicity. The laminar flame 
studies of IChamulak et al.] ( |200 8) predict a 10% reduction in 
the DDT density for a increase in the ^^Ne abundance of 0.02. 
From results found here, we estimate that this might corre- 
spond to a 3% decrease in NSE mass. This is about half as 
strong as the dependence from inherited neutron excess. 

Energy release - The amount of energy release per unit 
mass can affect both the global dynamics of the expansion 
and eventual disruption of the WD, as well as the more local 
dynamics related to the buoyancy that drives the acceleration 
of the burning front during the deflagration phase. The com- 
position is essential for both these effects in two ways. First, 
the gross nuclear energy available per unit mass is mostly sen- 
sitive to the abundance of '^C. Due to mixing during the smol- 
dering phase, the abundance over a broader region of the WD 
is involved in this case, as opposed to precisely at the cen- 
ter as is the case for the ignition density above. Thus, the 
products of both the helium core and shell burning are im- 
portant. Secondarily, the -^Ne abundance, by changing the 
balance of protons and neutrons, can influence the NSE state 
to which material will flash. More neutron rich material fa- 
vors m ore tightly bound m aterial and therefore releases more 
energy (|C alder et al.|2007| l. 

Flame speed - The propagation speed of the subsonic burn- 
ing front propagated by thermal diffusion, the flame, is sensi- 
tive to both the '^C and ^^Ne abundances due to their effects 
on the energy release and the speed of the early stages of the 
nuclear fusion (Timmes & Woosley| |1992{ Chamul ak et al.| 
|2007| l. This laminar propagation speed is likely to be much 
less important after the strong turbulence that results from 
the buoyancy instability develops. Turbulence has the abil- 
ity to add a nearly arbitrary amount of local area to the flame 
surface, therefore making the burning effectively independent 
of the laminar propagations speed (Niemeyer & Hillebrandtj 
1995| l. However, there is a period of time early in the defla- 



gration phase of the supernova where the interaction of the 
laminar flame speed with the lower strength turbulent veloc- 
ity field from the core convection will be important for setting 
the morphology of the burned region at the point when strong 
buoyancy takes over ( Zingale & Dursi,2007 ). In this work, we 
find an important sensitivity to the assumed outcome of these 
earliest stages. 

WD structure - The WD is supported by pressure of de- 
generate electrons. The neutron excess, and therefore the 
^^Ne abundance, sets the number of baryons worth of weight 
each electron must support. Thus a WD of the same mass 



with a higher neutron excess will be more compact because 
it will have fewer total electrons. (Conversely at the same 
central density, a star with more ^'Ne, and thus more neu- 
trons, will be slightly less massive.) This has a concomitant 
effect on the densit y distrib ution throughout the star This is a 
fairly small effect ( Hoflich et al.^1998) , but important to treat 
appropriately due to the marginally bound nature of a near- 
Chandrasekhar mass WD. 

Neutron excess - In addition to the indirect effects men- 
tioned in relation to energy release above, the neutron ex- 
cess, as set principally by the ^^Ne content, has a direct in- 
fluence over the final nucleosynthetic prod ucts, particularly 
the amount of ^^Ni. Timmes et al. (2003i showed that the 



decrease in the ^^Ni produced in the explosion, absent other 
factors, should be linearly proportional to the ^^Ne content 
and therefore the original metallicity of the stellar population. 
The distribution of ^^Ne will again be important, with that in 
the inner bulk of the star being important for the gross yield, 
and that in surface layers for opacity sources in those regions 
of the ejecta. 

In order to bring the scope of our study within tractable 
limits, it is necessary to either neglect or exclude some of 
these effects with assumptions. Most notably we would like 
to avoid for now the possibly quite complex dependence of 
ignition density on composition. This can naturally be done 
by considering only a single value of Xuq, and assuming that 
we are only studying the variation for progenitors that resulted 
from the same formation history. Note that since metallicity 
changes main sequence lifetimes, for example, this simplifi- 
cation may be a rather unnatural one in contrast to compar- 
ing SN la at the same stellar population age (|H6flich et al. 
|1998| l. Our second simplifying assumption will be to neglect 
the dependence of the DDT density on composition. Explor- 
ing the effect of DDT density will be undertaken in immedi- 
ate future work along with the compositional inhomogeneity 
with which it is intertwined. These two sets of assumptions 
leave us to study the dependence arising from how changes in 
^^Ne abundance modify the energy release, flame speed, WD 
structure and neutron excess. Because [Timmes et al.|p003] l 
have provided an excellent discussion of the direct influence 
of neutron excess, we focus mainly on how these other effects 
may modify the robust conclusion they reached. 

3. IMPROVEMENT AND EXTENSION OF NUMERICAL MODEL 

The essential components of the numerical code used in this 
work were presented in Townsley et al. (j2007 ) and work ref- 



erenced therein, but we will give a brief overview. The over- 
arching code is formed by the FLASH Eulerian compressible 
hydrodynamics code ( [Fryxell et al.|2000 Calder et al.|2002| l 
with modules added for the nuclear burning which occurs in 
SNe la. FLASH uses a high-order shock-capturing compress- 
ible hydrodynamics method, the piecewise parabolic method 
(PPM, [Colella & Woodward||1984|l, adapted to tre at a gen- 
eral equation of state (EOS, [Colella & Glaz|[T985l ). We use 
flash's tabulated fully ionized electron-ion plasma EOS 
([Timmes & Swestyl [200n} [Fryxell et al.|[2000| . FLASH ap- 
plies this hydrodynamics method on an adaptive ly refined, 
tree-structured, non-moving Eulerian grid. We make exten- 
sive usage of this adaptive mesh refinement (AMR) capabil- 
ity, using different resolutions for burning fronts (4 km), the 
initial hydrostatic star (16 km), and the region of negligible 
density initially outside the star (as coarse as thousands of 
km). Extensive detail on our refinement scheme and tests of 



resolution are given in section 3.3 below. 
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The nuclear burning processes, beginning with carbon fu- 
sion, and extending to electron capture in material in NSE 
are implemented with a nuclear energetics model described 
by Townsley et al. (2007), and here in section 3.1 and cali- 
brated to reproduce the features of nuclear processes t hat oc- 
cur in SN la as calculated using hundreds of nuclides (Calder 



et al. 2007 1. This burning model is used for both subsonic 
(deflagration) and supersonic (detonation) burning fronts. Be- 
cause the flame physics that governs the propagation speed 
of the subsonic burning front is unresolved at 4 km (Chamu- 
lak et al.|2008 1, we propagate this front by coupling our nu- 
clear energetics to an artificial, resolved reaction front given 
by the advection-reaction-diffusion (ARD or ADR) equation 
( |Khokhlov|1995[|Vladimirova et al.|2006| l. A variety of spe- 
cial measures are necessary to ensure this coupling is acous- 
tically quiet and stable, and therefore appropriate for simu- 
lating the buoyant instability of the burning front (Townsley 
|et al.|[2007] l. Detonation fronts are handled somewhat natu- 
rally by the shock-capturin g features of the hydrodynamics 
scheme ( Meakin et al.|2008) l. The only common components 
between our code and that of| Plewa|(200 7 ) are those publicly 
available as components of FLASH, which excludes all com- 
ponents treating the nuclear burning; differences are discussed 
in [Townsley et al. ( 2007 1. 

In this section we discuss various changes made to im- 
prove the overall consistency of the numerical modeling as 
well as the extensions necessary for the burning model and 
flame speed treatments to include an arbitrary ^^Ne content. 
Changes to the burning model, flame speed and refinement are 
discussed in turn, with separate subsections for improvements 
and extensions. The improvements include a backwards- 
differenced time integration method for the nuclear energet- 
ics model, more well-defined methods for defining the den- 
sity used to compute the flame speed and the Atwood num- 
ber used to calculate the buoyancy-compensated ("turbulent") 
flame speed, as well as simplified and better tested mesh re- 
finement prescriptions. Extensions include modifications to 
the coupling of the RD front and nuclear energetics in order 
to better account for the interaction of the detonation front 
with the artificially thickened sub-sonic burning front, and in- 
clusion of the dependence of the laminar flame speed on com- 
position. 

3.1. Burning Model 

Although the model of nuclear energy release used here is 
functionally identical to that of Townsley et al. ( 2007| l, it is 
useful to restate it in a way that makes adjusting the abun- 
dance of the fuel easier and that is more amenable to back- 
wards differencing or sub-stepping in time. This will also pro- 
vide the context in which we can present our prescriptions for 
how the detonation interacts with the artificially broad flame 
front. 

3.1.1. Simplified Dynamical Equations 

The first step is to abstract the properties of the fuel and the 
ashes of carbon burning into adjustable parameters. The prop- 
erties of interest are the number of electrons per baryon, Y^, 
the number of fluid ions per baryon, Fion, and the average nu- 
clear binding energy per baryon, q. These can be constructed 
from mass fractions Xj per 



1 



(2) 
(3) 



where Z, is the nuclear charge. A, is the atomic mass num- 
ber (number of baryons), and Eh.i = {Zimp-Nim„-mi)c~ is the 
nuclear binding energy, where A^, = A, -Z, and nip, m„, and m,- 
are the rest masses of the proton, neutron and nucleus of nu- 
clide i. Note that because the quantities , Fion and q are "per 
baryon", it is more clear to think of the density on the grid as 
representing the baryon number density. The advantage being 
that, unlike mass, this is actually a conserved quantity, making 
it possible to calculate the overall change in rest mass energy 
in a well-defined way. Baryon number divided by Avogadro's 
number makes an extremely good approximation for the mass 
in grams, and this correspondence will be used so extensively 
that in situations where the difference is of no consequence, 
we will almost always refer to the mass density instead. In 
this interpretation the X, are actually baryon number fractions 
and also become conserved quantities in the absence of trans- 
formations. 

We begin by defining the properties of "pure" fuel or 
carbon-burning ashes as Ygj,Yionj,qf and Ye,a,Yion^a,qa re- 
spectively. As before, we also define three progress variables, 
which measure the progress of the burning through various 
stages. These are for processing of fuel to carbon-burning 
ashes, (l)aq, for processing of these ashes to Nuclear Statistical 
Quasi-Equilibrium (NSQE), and finally (p^j^ for relaxation of 
NSQE to full nuclear statistical equilibrium (NSE). Finally 
we will define the properties of the NSQEh-NSE material, 
SYe^qNi^Yion qfjjSqgN in such a way that the bulk properties are 

Ye = (l-(l> faWeJ + (0/fl " + 5Y,^qN (4) 

Yion = (1 -4>fa)Y\oaJ + {4>fa " 4>aq)Y\on.a + ^Yioa^qN (5) 
q={\-(t)fa)qf + {<i>fa-(t)aq)qa + 5qqN ■ (6) 

T he usage of, for examp le, 5qqN rather than qq^, as was done 
in [Townsley et al. ( 2007| l, is to avoid issues of how the evolu- 
tion of the NSE material should be treated when (pq^ is very 
small. 

As discussed in 'Townsley et al. (2007) it is possible to cal- 
culate a "final" NSE state by calculating the endpoint of an 
isobaric or isochoric burning based on the instantaneous local 
conditions. If we denote this state by "NSE" and also use the 
NSQE and NSE timescales parameterized based on the ex 
pected temperature of this final state, t^sqe and ti^se (Calder 
|et al.|2007^ , we can posit Lagrangian source terms. 

D^fa 



(1) 



Dt 

D4>aq 

Dt 
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Dt 

D(5qqN) 

Dt 
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Dt 
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Dt 



= max(0, 
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'cc 
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Here (^rd is the contribution from the reaction-diffusion 
mode l for the burning front propagation ( [Townsley et al.| 
[2007] l. 

0rd = kV20rd+-7;((/.rd) (13) 

T 

where (/)rd is the progress variable for the (A)RD front and 
the coefficients k and r are determined from the prescribed 
propagation speed of the front, S, and its desired width. 
Ricj)) = /(0-e())(l -(/)+ei)/4, with tunable parameters /, and 
< eo,ei ^ 1, provides a stable and acoustically quiet reac - 
tion front propagation as discussed in 'T ownsley et al.| ( |2007| l. 
Thermal fusion is included via the reaction term 



(f>CC = pXnc,f(l-(l)faf^Nf,{av)c^c(T) 



(14) 



where Na{o'v)c+c(T) is the carbon-carbon fusion rate from 
Caughlan & Fowler ( 1988), p is the local density, and Xncj is 



the carbon mass traction m pure fuel. While this rate is sub- 
ject to several important uncertainties as described in section|2] 
above, the features of the detonation are fairly insensitive to 
its precise value, depending mainly on the overall energy re- 
lease. Some additional prescriptions for treating the temper- 
ature used to calculate this rate are required in regions where 
the RD front is active, (/)rd > 0, which will be discussed be- 
low. In order to accurate ly propagate the de tonation, (j)cc is 
set to zero inside shocks (Fryxell et al. 1989). Finally, Fion,^ is 
a very rough estimate of the ion abundance for NSQE material 
given by 



Yion.q - max 



1 

28 



mm 



'ion,NSE~ 
4 56 



5^1 




Conservation of mass energy yields the evolution of the mass- 
specific energy due to these nuclear processes. 



De Dq 

— =Na — - 
Dt Dt 



4'qN[Ye.NSENA(mp + - m„)c^ + e^^NSs] , (16) 



where Dq/Dt is obtained using (j6| and the above source 
terms, nip, m„, and are the proton, neutron and electron 
masses respectively, c is the speed of light, and e,^,NSE is the 
energy loss rate to neutrinos due to weak processes in the pr e- 
dicted NSE state ( |Calder et al.|2007[[Seitenzahl et al.|2009l l. 

All the dynamical quantities in equations ([7|l-( 12 1 are mass- 
specific, or more appropriately baryon-specific, quantities. 
Their full hydrodynamic evolution is thus given, for a rep- 
resentative quantity, q, by 



djqp) 
dt 



= [-V-(qpv)]^ + 



Dq 



(17) 



where v is the velocity field. The portions indicated as 
H and B are treated in an operator split fashion, each 
acting consecutively, with H being the conservative hy- 
drodynamic operator implemented with Piecewise Parabolic 
Method (PPM) (Col ella & Woodward |19"84| and B the burn- 
ing scheme described above. Because the B operator also does 
not change p, equations (|7])-([T2]) reduce to simple time deriva- 
tives after operator splitting. For the application of the B op- 
erator, we further assume that the temperature used in the cal- 
culation of 0CC as well as the timescales and properties of the 
NSE final state are constant, so that equations (|7])-([T2]) can be 
backwards-differenced and solved algebraically. This is in- 
tended to increase the stability of the treatment, particularly 



when tnse is similar to the size of a hydrodynamic time step. 
While the assumption that T is constant is not a good one for 
the thermal term, this is likely to have fairly little impact as 
this term is largely only active in propagating the detonation, 
for which the burning length is unresolved. The NSE final 
state prediction should only vary slowly with time, even in 
burning fronts, since it depends on the combination e-NAq or 
e-PNAq/ p for the isochoric and isobaric predictions respec- 
tively. The former combination only varies on the hydrody- 
namic timescale and the timescale of weak processes, and the 
latter should vary fairly slowly due to the subsonic propaga- 
tion of the RD front. 

3.1.2. Detonation- Flame Interaction 

The RD front that we utilize to model the propagation of the 
subsonic burning front is several computational zones thick 
in order to enable its stable propagation. This presents sev- 
eral complications when it is desirable to also treat thermally 
activated burning that is not related to the flame front. At 
high densities, the physical flame is very thin, and as a result, 
there is very little partially burned material. Thus, in this den- 
sity regime, partially burned material on the grid represents 
regions where burned and unburned material are well sepa- 
rated, but the interface is not resolved. In addition, this means 
that the temperature of the zone is not representative of the 
temperature in the unburned material, which s hould be used 
in the calculation of (/)cc- In previous work ( Meakin_et al. 
2008 I, the thermally activated term, 0cc in equation Q was 
suppressed (set to zero) within the RD front, where 0rd ex- 
ceeded some small threshold. This is not too bad an approx- 
imation in the case of the gravitationally confined detonation 
(GCD) scenario because only a small portion of the star is 
burned via the deflagration, so that an even smaller portion is 
left in an unrealistic partially burned state. However, such a 
prescription would leave a large amount of partially burned 
material in a simulation in the DDT scenario. 

The principal intent of allowing a thermal contribution 
within the RD front is to allow the detonation to consume 
nearly all of the partially burned material as it encounters the 
RD front, as it would if encountering a thin flame. In order to 
accomplish this, two measures are applied. First, 0cc is sup- 
pressed where (p-^D > 10"'' unless ^/a-^rd > 0.1 within one 
flame width (4 zone widths) of a given cell. This prevents spu- 
rious thermal burning. Second, where (pcc is enabled within 
the RD front, an estimate of the temperature of the unburned 
material is used in place of the local temperature in the calcu- 
lation of {(Tv)c+c{T)- This temperature estimate is obtained by 
estimating the properties the local material would have in the 
absence of the RD-front based fuel consumption. The frac- 
tion burned due to thermal processes is Xh.th = ^/o-^rd- The 
properties in the absence of the RD front would then be 



th = (1 -Xbj)Yion.th +Xb 



where we estimate 



qb = 



Y\f\n h — 



q-(l-(j)fa)qf 



Ymn-i^-4'fa)YionJ 



(18) 
(19) 



(20) 
(21) 



Then the temperature estimate for the thermally activated 
burning, 7ih, is found by evaluating the equation of state for 
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the local density, Y^, Fion.th and mass-specific internal energy 
given by e± = e-{q-qth)- With a sufficiently accurate and ro- 
bust estimate of the temperature, the first condition, suppress- 
ing (j>cc in the RD front except in proximity to other burning, 
might prove unnecessary. This was not the case in tests we 
have run so far The above prescriptions appear sufficient for 
the purposes of this study. 

3.2. Flame Speed 

The improvements to the flame speed treatment are in- 
tended only to improve numerical consistency by obtaining 
a more well-defined and stable value for the front propagation 
speed. The improvements are to estimate the unbumed den- 
sity and to use a tabulation of the Atwood number instead of 
a crude estimate. In doing this, we also develop methodology 
for testing the accuracy of tabulations to ensure that the tabu- 
lation grid is sufficiently fine. Following this, we present de- 
tails of our extension of the computation of the laminar flame 
speed to account for dependence on fuel composition. This 
is one of the physical effects that our study of Ne systemat- 
ics is evaluating. The tabulation utilized for the laminar flame 
speed is also evaluated for accuracy. 

3.2.1. Consistency Improvements of Flame Speed Treatment 

To be consistent with the burning model, the input flame 
front speed is a function of the local chemical composition of 
fuel and the density of that fuel. We define the flame front 
speed to be the speed of the carbon burning front with respect 
to the fuel. We parameterize the chemical composition of the 
fuel by '^C and ^^Ne mass fractions, with the remaining ma- 
terial being "'O. In order to prevent the flame from being torn 
apart by turbulence induced by the Rayleigh-Taylor (RT) in- 
stability, we introduce a minimum flame speed. 



5min = 0.5VAgmA, (22) 

where A is the Atwood number for carbon burning, g is the 
gravitational acceleration, A is the grid resolution, an d ot is a 
caUbrated parameter determined to be m ~ 0.04-0.06 ( Towns- 
ley et al.|2008p . In this work, we use m = 0.04. The Atwood 
number and gravitational acceleration describe the strength of 
the RT-induced turbulence. The resulting speed used to set the 
propagation of the ADR front is 



S — max (5niin , ^lam) , 



(23) 



where 5iam is the speed of the laminar flame front. S is then 
rolled off to zero for densities below 10^ g cm~^. In reality, the 
RT-induced turbulence wrinkles the flame, increasing the sur- 
face area of the burning front and accordingly increasing the 
burning rate. The minimum flame speed for our model thus 
serves to compensate for this boosted burning rate jKhokhlovl 
1995| l. By construction, this prescription does not properly ac- 



count for the interaction of isolated turbulence with the flame 
surface. Thus while the integrity of a given plume is main- 
tained, the interaction of turbulence created by one plume 
with others is not necessarily captured accurately. We per- 
form resolution studies and vary the parameter m to check for 
precisely these issues (Section \53\ and find that they do not 
appear to adversely affect the principal metric measured in 
this study, the expansion prior to DDT. Inclusion of more di- 
rect models of flame-turbulence interaction ( jColin et al.|2000t 
[Schmidt et aLp006| l is planned for future work. 

Both 5iain and A are best characterized by a dependency 
on the composition and density of the unburned fuel. This 



presents a difficulty when S must be constructed in a partially 
burned cell, where the density and composition has changed 
from that of the unburned material due to burning and energy 
release. For the composition, we simply store separate mass 
scalars that represent the initial composition: the mass frac- 
tions of material that was initially in the form of '^C and ^^Ne. 
While the density varies across the subsonic burning front, the 
pressure is approximately constant. We therefore use the lo- 
cal pressure, P, to estimate the density of the unburned state. 
This is greatly simplified by the fact that the unburned ma- 
terial in the WD is at a high degree of degeneracy, so that 
the temperature of the unburned state is not too important and 
simple forms for the pressure-density relation of a degenerate 
gas can be used (Ha nsen & Kawaler|I994 i. Our estimate for 
the unburned density is 



Pu = 0.92 X 




1.243 X 101 



3/2 . p .6/5 



1.004 X 101 



(24) 

where P is the pressure in cgs units (erg cm ^), is the elec- 
tron fraction, and 0.92 is adjusted to provide a good fit to the 
pressure-density relation in the initial star. The difference be- 
tween this estimated density based on the local pressure and 
the actual density is no more than 6% in the initial star for 
p > 2 X 10^ g cm"^ and varies smoothly from an overesti- 
mate at low pressures to an underestimate at high pressures. 

The Atwood number used in the calculation of the mini- 
mum flame speed is determined based on the local fuel com- 
position and the resulting energy release from carbon burning. 
Given these parameters, we estimate the density of the ash 
using the Rankine-Hugoniot jump condition across the flame 
front ( [Vladimirova et al.|2006 1. For computational efficiency, 
we calculate the Atwood number at the beginning of the simu- 
lation for twenty-one equidistant log densities between 6 and 
9.6 and ten equidistant carbon mass fractions between 0.3 
and 1.0 using a representative neon mass fraction appropri- 
ate for the simulation. Because the Atwood number changes 
less than 0.01% over the relevant range of X22, we reduce the 
dimensionality of the interpolation and memory requirements 
by introducing a characteristic neon mass fraction, which for 
this study is equal to the global neon mass fraction. During a 
simulation, the code performs a bi-linear interpolation of the 
Atwood number given the local initial mass fraction of '^C 
and the estimated unburned density from Eq. p4| . 

In order to test the accuracy of this interpolation procedure, 
we estimate the uncertainty for the Atwood number through- 
out the table. Our method of uncertainty estimation solves for 
the second-order term in a polynomial interpolation of the ta- 
ble in each direction. Because the Atwood number is a smooth 
and slowly varying function of X12 and logjQp as shown in 
Figure[T] we can assume that higher order terms in the Taylor 
expansion are negligible. In this case, the second-order term 
serves as a correction to the linear interpolation in a particu- 
lar direction. We use this correction to define the uncertainty 
estimate 



Ri(Xi2Aogiop) 




(25) 



where the Atwood number, A = A(Xi2,logjQp), is calcu- 
lated using a 1st and 2nd order interpolation. The second- 
order interpolation, A^"*^, is performed in only one variable- 
dimension, ^= {Xi2,logiQp}. By subtracting the first-order 
from the second-order interpolation, we are left with a term 
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Fig. 1. — The Atwood number is shown as a function of log density for 
several carbon mass fractions. This displays the sufficient smoothness of the 
Atwood number as a function of log density and carbon mass fraction such 
that our method of uncertainty estimation is valid. 

proportional to the second partial derivative with respect to 
^. This allows us to assess the curvature of our table in each 
direction and whether we have enough points in our table in 
each dimension such that linearly interpolating the data pro- 
vides an accurate estimate of the Atwood number. For the 
Atwood table, we are able to calculate an exact value to com- 
pare against the interpolated estimate. However, for the flame 
speed table discussed later, we are not able to calculate an 
exact value to test against. We use the Atwood table to ver- 
ify our method of uncertainty estimation used later to analyze 
our flame speed table. The comparison of the uncertainty esti- 
mate, to the actual error, a, in Table [l] verifies our method 
of uncertainty estimation. 

The results of performing these tests at the grid points and 
midpoints of the Atwood table are shown in Table[T] We show 
the minimum, maximum, absolute average, and average un- 
certainty estimate, Rmax-.Rmin, (1^1)^ in each direction of 
the Atwood table. In columns 1 and 2, we evaluate the uncer- 
tainty estimate at a grid-edge (a midpoint along ^ holding the 
other variable at tabulated values). In column 3, we evaluate 
the uncertainty estimate at a cell-centered point (where both 
variables are at midpoints). In this case, we sum the uncer- 
tainty estimates in each direction, R^^^ = Rxi2~^^iogp- '^^ 
not perform a bi-quadratic interpolation. While the cross-term 
is potentially important in this case, we do not include it in the 
uncertainty estimate of a cell-centered point due to the depen- 
dence on the order of interpolation in multi-dimensional poly- 
nomial interpolations. In columns 4 and 5, we calculate the 
actual error, a, for any X12 and X12 = 0.5, respectively. In this 
study, we use models withXi2 = 0.5. The error is calculated at 
all midpoints and grid points (the same values used for the un- 
certainty estimates). Due to the way we initialize the Atwood 
table, the test points are equivalent to forty-two equidistant 
log-densities between 6.0 and 9.6 and twenty equidistant car- 
bon mass fractions between 0.3 and 1.0 with X22 = 0.01. 

Because the table has positive curvature in X12 and negative 
curvature in logj,,/?, we expect to find the maximum magni- 
tude in the uncertainty estimate on a grid-edge as opposed to 
a cell-centered point. Our method of uncertainty estimation 
provides a maximum expected uncertainty of -10.99% which 
is verified by the true maximum error of -9.99%. Both of these 
values were evaluated at X12 = 0.339 and logjQ p = 9 All. This 
shows that our method of uncertainty estimation works for the 



TABLE 1 

Error estimates for tabulating and linearly interpolating 
THE Atwood number 



€ 


Xn (%) 


logioP(%) 


Xl2+logioP(%) 


a (%) 


O"C12=0.5 {"/") 




-0.98 


3.30 


0.00 


0.03 


-0.15 


^max 


-10.99 


6.88 


-5.19 


-9.99 


-7.60 


(l«l> 


4.82 


5.52 


1.93 


4.11 


3.69 


{R) 


-4.82 


5.52 


0.82 


0.95 


-3.10 



Atwood table and can be applied to the flame speed table that 
will be discussed next. For the purposes of this study, the rel- 



evant maximum error is cr, 



ci 2=0.5 



: 7.6% obtained from com- 
paring the interpolation against calculated Atwood numbers 
at midpoint log-densities all at X12 = 0.5. 

3.2.2. Composition Dependence of the Laminar Flame Speed 

For the laminar flame speed, we give preference to values 
calculated by Chamulak et al.| (2007) using a 430-nuclide re- 
action network for a variety of initial carbon and neon mass 
fractions and a range of densities. Similarly to the Atwood 
numbers, at runtime the code performs a linear interpola- 
tion to obtain the laminar flame speed from a table of previ- 



ously calculated results. The method used by Chamulak et al. 



2007[ l is not well suited to solve for the laminar flame speed 



at low density; therefore, we use the results from Timmes &| 
|Woosley| ( |1992| l to supplement this table. Because the flame 
speeds for each case were calculated using different initial 
carbon mass fractions, we merged the data by linearly interpo- 
lating Tim mes & Woosley| ( |1992| l values onto the Chamulak] 
iCt al., ( ,2007) 1 grid. The results of this merger are more clearly 
shown in Figure [2] 

Some discussion of the accuracy of this method of ob- 
taining flarn e speeds at low densities is warranted. |Timmes| 



& Woosley 



( 1992| l do not track ^^Ne dependence with their 
method of determining the laminar flame speed. For identical 
points in parameter space, the two methods produce laminar 
flame speeds that differ on average by ^ 30%. This differ- 
ence is of the same order as the effect of adding ^^Ne where 
we see a ^ 30-60% speed-up depending on density. For the 
models considered in this study, low densities occur near the 
surface of the white dwarf star In these regions, the input 
flame speed (23 1 is dominated by the Rayleigh-Taylor driven 
turbulence such that S'niin > •S'lam- In fact, this transition oc- 
curs at /9 sa 2.5 x 10^ g cm"-' in the initial star. Therefore, 
this tabular method of estimating the laminar flame speed is 
sufficient at low densities for this study. 

The tri-linear interpolation occurs within the three- 
dimensional parameter space of carbon and neon mass frac- 
tion and log-density to calculate the laminar flame speed. 
We cannot easily evaluate the uncertainty in our interpola- 
tion method by comparing with a direct calculation of the 
laminar flame speed due to the computational cost. There- 
fore, we apply our method of uncertainty estimation discussed 
for the Atwood number to the laminar flame speed table, 
Eq. (25 1. We calculate an uncertainty estimate, for our 
methoa of interpolation at the quarter-points and grid-points 
for each parameter in our table. We limit our analysis to den- 
sities above 2.5 x 10*^ g cm"-*. The laminar flame speed 
becomes unimportant below 2.5 x 10** g cm"-' because the 
buoyancy-compensated flame speed takes over in Equation 
( 23 I at roughly this density. The results of these calculations 
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Fig . 2 . — Shown are l aminar flame speeds as a fun ction of log density from 
IChamulak et al. i'2007), [Timmes & Woosley|rr992) , and our resulting flame 
speeds from merging these datasets all at zero ^Jemass fraction. Each panel 
compares these datasets at different carbon mass fractions, X|2. *(Note that 
ITImmes & Woosley ( 1992 1 performed their study at Xyi — 0.2. We expect 
slightly higher laminar flame speeds atXi2 = 0.3.) 

are given in Table |2] and Table |3] Table [3] shows results for 
X\2 = 0.5, which is more relevant for this study, while Table|2] 
shows the behavior of our flame speed table in general. As 
was the case for the Atwood number, the uncertainty estimate 
along a single variable dimension, is evaluated along a 
grid-edge (at the quarter-points of ^ holding the other vari- 
ables at tabulated values). This is the case for columns 1-3. 
In column 4, we evaluate Xi2-log[Q/o cell-centered points (at 
the quarter-points of X12 and logjoP holding X22 at tabulated 
values). To evaluate the uncertainty estimate on cell-centered 
points, we sum the individual uncertainty estimates in each 
direction. In column 5, we calculate the uncertainty estimate 
at the quarter-points in all directions such that X12, X22, and 
logjop are not at tabulated values. 

For the flame speed table, both the X12 and logjgp variable 
dimensions have negative curvature, while X22 has mostly 
positive curvature. Because the uncertainties along each grid- 
edge are comparable, we expect the maximum uncertainty in 



TABLE 2 

Error estimates for tabulating and linearly interpolating 

THE LAMINAR FLAME SPEED (ABOVE log,o p = 8.4) 





Xi2 (%) 


X22 (%) 


logioP(%) 


Xi2+log,oP(%) 


X,2+X22+log,oP(%) 


^min 


0.00 


0.00 


0.24 


0.45 


0.39 


^max 


8.74 


-10.96 


5.81 


8.63 


8.27 


m 


1.14 


0.25 


1.76 


2.73 


2.68 


(«> 


1.01 


-0.22 


1.76 


2.73 


2.68 



TABLE 3 

Error estimates for tabulating and linearly interpolating 

THE LAMINAR FLAME SPEED (ABOVE logjo p = 8.4 AND X[2 = 0.5) 





X22 (%) 


logioP(%) 


Z22+log,oP(%) 


^min 


0.00 


0.70 


0.72 


^max 


0.62 


3.35 


3.34 


(l«l> 


0.06 


1.82 


1.83 


(R) 


0.06 


1.82 


1.83 



our table to occur either along a grid-edge or on a Xi2-log[g p 
cell-centered point. We do not expect a maximum uncertainty 
estimate involving X22 and any of the other two variables due 
to opposing curvatures. 

We determined that the laminar flame speed table has suffi- 
cient resolution at densities above logjQ/9 = 8.4 with the mag- 
nitude of estimated uncertainties < 10% as shown in Table |2] 
The estimated uncertainties relevant for this study at X[2 = 0.5 
are < 3% as shown in Table [5] 

3.3. Mesh Refinement 

The goals of the design of our refinement scheme are to be 
as simple as possible while both capturing interesting or phys- 
ically important features and doing so with good efficiency. 
These three goals are somewhat at odds, and therefore pro- 
vide a wide latitude for choosing refinement prescriptions. We 
proceed by defining three regions of the physical domain: 

1 . fluff: regions with p < pfiuff 

2. star: non-fluff regions, p > pfluff 

3. energy generation: 

regions with e^uc > feg or 0rd > (j^Ro.eg 

where (/)rd is the progress variable in the reaction-diffusion 
model of the flame front. These are indicated with a "/", "*" 
or "eg" subscript respectively. We assign to each of these a 
consecutively increasing maximum refinement level. For sim- 
plicity, we will here use the minimum cell size rather than the 
refinement level. A/ > A* > Agg. 

A "fluff" region outside the star is necessary because the 
hydrodynamics method in FLASH has no explicit mechanism 
for treating empty (zero-density) cells or free surfaces. To 
ameliorate this, would-be empty cells are filled with material 
of extremely low density which will not effect the dynam- 
ics of th e more dense materi al of interest in the simulation 
(see also Zingale et a l.|2002 i. Here we set it to 10"^ g cm~^. 
Because this material is of no physical interest and has neg- 
ligible contribution to the dynamics of other material, A y is 
taken as large as possible, generally being the total domain 
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size divided by the block size. Recall that the smallest inde- 
pendently refinable unit in PARAMESH is a "block", which 
in all our simulations is a 16 x 16 cell region. 

The finest resolution Aeg will be treated as the resolution 
of the simulation, though in fact before the burning spreads, 
most of the star is only refinable to A*. Some care must 
be taken in defining the thresholds such that regions near the 
thresholds do not cycle between refinement levels. Cycling is 
avoided by checking both parent and child blocks (i.e. blocks 
at the finest and next to finest refinement levels) reconciling 
the results and using some hysteresis in defining the bound- 
aries. For the simulations presented here pfiuff = 10^ g cm"^ 
throughout the simulation. If a simulation is run far enough 
into expansion, this restriction would need to be relaxed to 
lower density so that stellar material remains refined as it 
expands. The energy generating thresholds are placed at 

= 0.02 s"' for the deflagration 



eeg = 10' 



erg g ' s ' and 



and detonation phases, but are relaxed to 10 erg g s" and 
infinity, respectively, for calculations of the expansion phase. 

Beyond the definition of these regions of maximum re- 
finement, refinement is triggered by gradients in the physical 
quantities. Here we use the density and the first progress vari- 
able in the burning model, (f)f„, to trigger refinement. This 
choice guarantees that all burning fronts are refined to the de- 
gree allowable. In order to detect gradients, we utilize the 
built-in tests included with PARAMESH and FLASH. These 
measure a ratio between the second and first derivatives of the 
fields being checked. The reader is referred to the source code 
for the exact generalization to multiple dimensions. A thresh- 
old for triggering refinement is then defined for each quantity 
checked, the default being 0.8. In order for the initial star to be 
fully refined, it was necessary to drop this threshold to 0.1 for 
the test of the density field. With this threshold, de-refinement 
set in during expansion in one-dimensional simulations, but 
not appreciably so in two dimensions. 

In order to verify that sufficient resolution was being uti- 
lized and that the adaptive refinement was not adversely af- 
fecting the accuracy of the solution, a one-dimensional con- 
vergence test was performed. Convergence is expected in one 
dimension because instabilities are suppressed. The absence 
of the instabilities that accelerate the burning in a multidi- 
mensional simulation necessitates a significant artificial en- 
hancement of the burning rate in order to unbind the star. Be- 
cause this is only a numerical test, we simply increased the 
m parameter (Equation 1 22 1) until an explosion was obtained. 
Then the combination inK, where A is the finest resolution, 
generally Aeg, is held fixed as the resolution is varied. We 
used inA = 32 km. We are not interested in the order of con- 
vergence here, only that reasonable convergence is obtained. 
Therefore, we simply compare the density distribution of the 
outgoing ejecta to a uniformly refined case with additional 
levels of refinement, Aeg = A* = 1 km. This is the solid line 
in Figure[3] which shows density distributions for simulations 
at various resolutions at a time of 5.6 seconds after ignition. 
The fluff was not refined in any of these simulations, but a 
test was performed with the entire domain refined to confirm 
that this had no effect in the density range shown here. The 
same initial WD was used in all cases, which was created on 
an 8 km grid and mapped onto the new grid by averaging den- 
sity and temperature without a reconstruction (or equivalently 
a piecewise constant one). 

We found that a resolution of Aeo = 4 km and A* = 16 km 
was necessary to satisfactorily match the reference solution. 




radius (10 km) 

Fig. 3. — Density distribution of outgoing ejecta in radius for one- 
dimensional test of resolution. Convergence is obtained when energy- 
generating regions are refined to 4 km and the rest of the star is initially 
refined to 16 km (dashed line). The reference case in which the entire star is 
uniformly refined to 1 km is shown by the solid line, while cases in which 
each of the refinement regions are relaxed by one level from the coarsest con- 
verged values ai'e shown by the dot-dashed and dot-dot-dashed lines. 

Also it was necessary for the initial star to be fully refined at 
A*, and the trigger threshold for the density gradient test was 
adjusted to achieve this as described above. Subsequent de- 
refinement as the star expanded did not appear to have any 
adverse effect on the outcome of the simulation. The de- 
refinement threshold was set to 0.0375. The case with these 
parameters matches the reference case extremely well down 
to a density of about 3 x lO-' g cm"^. Most likely the initial 
resolution is insufficient to resolve the hydrostatic equilibrium 
in these low density layers. For comparison, we show in Fig- 
ure[3]the cases with a factor of 2 too little resolution separately 
for the energy generation region and the star. We find the out- 
come is most sensitive to the refinement level of the star. If 
this is too low, the hydrodynamics of the stellar expansion are 
not properly captured. 

This convergence study also provides an important check 
on the model of nuclear energy release, particularly demon- 
strating that sub-stepping between hydrodynamic steps does 
not appear necessary for such a simple reaction model, as it 
likely would for even a highly reduced nuclear network. It is 
possible that lower resolutions of the energy generation region 
could be made viable by introducing a sub-stepping mecha- 
nism. 

4. FRAMEWORK FOR EVALUATING SYSTEMATIC DEPENDENCIES 

As a starting point for our study of systematic effects, we 
must consider a model of the supernova explosion that repro- 
duces many observed characteristics. In light of its success 
in one-dimensional work (e.g. |Hoflich & Khokhlov||1996"l ), 
we take the deflagration-detonation transition as this starting 
point. The defining feature of this scenario is that the flame 
front will undergo a transition to detonation when turbulent 
mixing on the scale of the flame front becomes more vig- 
orous compared to the propagation speed of the flame front 
(|Niemeyer (fe Woosley||1997t [Khokhlov et al.||1997|l. One- 
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dimensional work such as |Hoflich & Khokhlov| ( |1996| l uti- 
lized the density at which the DDT occurs, pom, as a non- 
unique parameter. Variation of pddt then led to the observed 
variety of SN la outcomes. Instead of this, we make the as- 
sumption that the conditions that lead to the DDT, while de- 
pendent on local properties (e.g., composition and turbulence 
strength) are otherwise unique, though not currently known 
with precision. Variety in the outcome for the same initial 
WD relies on the ignition configuration being non-unique due 
to the turbulence present in the convective core at the time of 
deflagration ignition. In this section we construct the essen- 
tial component of our framework for evaluating systematic 
dependencies: a theoretical population of SNe la obtained by 
sampling a defined ensemble of ignition conditions. Some ba- 
sic characteristics of this population are discussed. The next 
section, [5] describes a typical outcome in detail along with 
additional technical details of the implementation of the ex- 
plosion. 

The ability of isolated plumes to rise to low densities from 
a central ignition before the star expands significantly, as 
demonstrated by a number of studies (ICalder et al. 1 120031 
Hewa et al. 2004; Calderet al. 2004; Livne et al. 2005; Plewa 
2007, ,Townsley et al.)|2007; |Ropke et al.||2007; ,Jordan et al.| 
2008| , presents a significant challenge to producing a realis- 
tic explosion with a DDT. If the transition to detonation takes 
place as early as these simulations imply, only a very thin 
layer of Si-group elements, the hallmark of a la spectrum, are 
produced. However, |Ropke et al. ( 2006, ) found that a multi- 
spot deflagration ignition led to a much more symmetric de- 
flagration. Building upon this, somewhat like other authors 
( |Ropke & Niemeyer 2007 ; Bravo & Garcfa-Senz 2008 ), we 
here are able to obtain realistic DDT explosions by consid- 
ering deflagration ignition conditions chosen to have no low- 
order power 

4.1. Constructing a Theoretical Population of SN la 

In order to evaluate systematic dependencies in the SN la 
population, we require a theoretical sample of supernovae that 
mimic the properties of supernovae as a population. SN la 
possess an intrinsically large scatter in the ^^Ni mass synthe- 
sized during the explosion, ranging typically between about 
0.3 and O.SM©. (See [Howell et al.|[2509bj for one example 



flame surface is defined to be 



r(9) = 150 km+30 km V AeY?(9) 



(26) 



sample distribution.) With the DDT scenario, we have found 
that this degree of scatter can be obtained by introducing a 
certain degree of randomness into the initial condition of the 
deflagration phase. We now motivate and describe this initial 
condition and then describe the samples that it produces. 

Our initial condition is motivated as a possible abstraction 
from prospective three-dimensional studies of the very early 
deflagration phase. During the first ~ 0.1 s of the deflagra- 
tion, the flame surface will be heavily influenced by the pre- 
existing convection field in the core. In order to develop a 
randomized sample of a variety of ignition conditions, we 
choose to parameterize the possible outcomes of this early 
evolution as harmonic structure in a flame surface when it 
reaches a modest distance from the core. We place the posi- 
tion of the initial flame surface before perturbation at a radius 
of ro = 150 km. To this base radial position, perturbations with 
definite harmonic content are added in the form of spherical 
harmonics. In the case of axisymmetry, only m = spheri- 
cal harmonics can be included, reducing to Legendre polyno- 
mials. However, this technique should extend in a straight- 
forward way to three dimensions. The initial position of the 



The normalization convention of [Jackson] ( | 1 99 9^ is used. The 
amplitudes are chosen from a normal distribution. The 
maximum harmonic content £n,ax is chosen so that the per- 
turbations in the flame surface are resolved. For the 4 km 
resolution simulation performed here, we choose ^^ax = 16. 
As demonstrated by a number of studies of locaUzed, off- 
center ignition in the absence of a convection field (|Calder[ 
et al. 2003 ; Plewa et al. 2004; Calder et al. 2004; Plewa 2007f 
,Ropke et al..2007,.Townsley et al.^2007,, Jordan et al..2008p , 
low-order perturbations on such a flame surface would lead to 
an early DDT with too little expansion of the WD to give re- 
alistic yields. To prevent this effect, low order modes are left 
out by choice of ^min. Here we have used £,^in = 12, which ap- 
pears to give a reasonable sample (see below). An extensive 
study of the sample that arises from different ^min values was 
not undertaken, so the sensitivity of the sample to this choice 
is uncertain. One additional restriction was imposed due to 
the axisymmetry. The perturbation, the second additive term 
in equation ( [26| , is restricted to be negative on the symmetry 
axis, 9 = 0, IT. This suppresses the formation of slender "jets" 
of burned material, which are likely artifacts of the axis sin- 
gularity. Such jets are much smaller than the rising bubbles 
observed by Townsley et al. (20071 and similar simulations, 
and bear more resemblance to the "tails" seen in those cases. 



Within the framework defined by equation (26 1, besides 
imin and imax, it in necessary to specify the implementation 
of the random choices of the Ag.'^ We will refer to a set of 
these as a realization of the initial condition. Rather than 
choosing completely unrelated random seed for each realiza- 
tion of the initial condition, a more well-defined and repro- 
ducible sample is obtained by drawing random numbers for 
each consecutive realization from a single random number 
stream. Thus the entire sample is represented by the initial 
seed of the random number stream, along with algorithmic 
details. Also this gives a definite ordering to the realizations, 
arising from the order of the random number stream. We do 
not require a large number of random numbers, and therefore 
opt for a simple linear congruential generator (LCG) pseudo- 
random number generator with a 31 -bit seed/output value. 
The ending seed from one realization is simply used as the 
starting seed for the next. Candidate realizations that do not 
satisfy the property of having a negative perturbation on the 
symmetry axis are dropped and another is generated. These 
dropped candidates are not included in the numbering of the 
realizations used in the n ext section. We use the LCG dis- 
cussed in section 16.1.3 of [Newman & Barkema[ ( [T999] l. The 
initial seed for the set of realizations presented here was ob- 
tained from the standard Linux kernel random number source 
(/dev/random) and is 1866936915. 

These sets of initial conditions provide a framework within 
which we may study a wide variety of possible systematic ef- 
fects. This includes both physical systematics such as those 
explored in this study, or those arising from physical uncer- 
tainties in the numerical model. As study of the central ig- 
nition mechanism for SNe la advances, with improved flame 
models and DDT conditions, for example, it may be neces- 

' Our implementation of the process described liere is available from 
|http : //variable ■ as ■ arizona ■ edu/code| 
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sary to adjust the overall parameters of this framework, no- 
tably ijnin, in order to maintain a realistic sample. Also, as 
discussed previously, these choices can be compared to simu- 
lations of the early deflagration phase in order to both inform 
future choices and retrospectively understand the context of 
previously performed studies. One advantage of such a con- 
trolled initial condition is that it provides a slightly stronger 
probe of systematic effects than would generally be available 
observationally. This arises from the fact that it is possible 
to compare the same ignition sample rather than independent 
samples, though the latter can also be performed if desirable. 

4.2. The Theoretical Population 

The above development provides a clear definition of a pop- 
ulation of ignition conditions from which we may draw a sam- 
ple. Lacking further information about the nature of the har- 
monic content of the initial condition, we will assume cases 
drawn from this sample will have equal weight (likelihood). 
The population of supernovae that results is a purely theoret- 
ical construct. Any observed population will have a variety 
of progenitors that will have a distribution of intrinsic proper- 
ties (composition structure, accretion history, etc). Therefore, 
while this population of ignition conditions will be very useful 
for studying systematic effects, caution should be exercised 
when drawing conclusions about observational ramifications. 
Future studies will address a more observationally motivated 
sample. 

In order to understand the diversity in our sample, and the 
qualitative changes arising from changes in ^^Ne fraction, the 
first 5 realizations from the sample sequence were run through 
the end of the detonation phase. The basic outcome proper- 
ties of these ten cases are listed in table |4] The DDT time, 
toDT, is defined as the time at which any part of the flame 
surface first passes a density of 10' g cm""*. This is the time 
at which the first detonation front is launched. The mass at 
high density at the DDT time, in this case that above a den- 
sity of 2 X 10' g cm"-', will be discussed more below. The 
NSE mass, Mnse, is defined as J (f)q„pdV integrated over the 
star. The ^^Ni mass is determined by using eq. (29 1 below to 
estimate the local mass fraction from the Yg and again integrat- 
ing over the whole ejecta. The Si-group and 0-Si masses are 
defined as J {(t)ag-<Pqn)pdV and J (i^a-4'atj)pdV respectively 
and are discussed more in section [Sj where the total nuclear 
energy released A/iiest mass is also defined. All of these yields 
are determined after the detonation phase is complete, and the 
gross production of burning products is complete (see section 
|5]). The amplitudes of the additive components that make up 
the initial conditions are given in Table[5] though this does not 
bear out any particular pattern in the results under discussion. 

Even this small sample from the theoretical population pro- 
duces a diverse set of supernovae. The estimated yield of 
""^Ni spans a range between 0.45 and 0.8Mq. The average 
is in the upper portion of this range. For Xiij^^ = the aver- 
age is O.7O±O.O5M0, with a sample standard deviation (a) of 
O.IIMq. ForX22Ne = 0.02 the average is O.62±O.O5M0, also 
with a sample standard deviation of 0.1 IMq. The decrease of 
the ^^Ni yield with increasing ^-Ne content is also directly re- 
flected in the fraction of the NSE material that becomes ^''Ni. 
This fraction is 0.86 with a = 0.01 for Xn^^ = 0, but is 0.80 
with a = 0.02 for X22^g = 0.02. This is a direct result of the 
differences in initial neutron fraction due to the presence of 
22Ne. 

The main determining factor in the gross amounts of prod- 
ucts synthesized in the explosions is the degree of expansion 
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Fig. 4. — Mass of all Fe-group (NSE, open symbols) elements and esti- 
mated mass of ''^Ni (solid symbols) produced in the explosion plotted against 
the mass of material above a density of 2 X lO' g cm"^ when the first deto- 
nation begins, an indicator of the degree of expansion of the star prior to the 
DDT. Two initial abundances of -^Ne are included, ^22^^ = (squares) and 
-f22Ne = 0.02 (circles). Both of these yields are well-con'elated with the mass 
at high density at the DDT transition, making the latter a useful intermediate 
indicator of the explosion outcome. 

when the detonation ignites. This dependence has been dis- 



cussed by previous author s (e.g. Ropke & Niemeyer 2007 
Townsley et aLf2 007; Me akin et al |2008] l 
such a relation will differ among different ( 



but the details of 
delayed detonation 
scenarios. It is useful to characterize the degree of expansion 
by the mass at high density (above some density threshold), 
which forms an indicator of how much material will be pro- 
cessed to NSE. A part of this will become ^''Ni and deter- 
mine the overall brightness of the supernova. Several density 
thresholds for defining the mass at high density were consid- 
ered. The mass above a density of 2 x 10' g cm"^, M(pj > 2), 
appears the most appropriate. This conjecture is demonstrated 
by the open symbols in figure |4j which lie near a 1-1 relation 
(dashed line). 

The variation in expansion at the DDT appears to arise from 
a competition between the rise of the highest plume and the 
expansion of the star in response to energy input. Both these 
processes have inherently similar timescales, the star's dy- 
namical time. Figure [5] compares the state of the star ap- 
proximately 0. 1 seconds after the launch of the first detona- 
tion for the 5 realizations at two ^^Ne fractions run through 
the detonation phase. The top row displays the Xiiyi,. = and 
the bottom Xnj^^ = 0.02, with each column being a separate 
realization of the ignition condition. The coloring indicate 
the nucleosynthetic yield, with black indicating NSE material, 
green Si-group material, and red 0-Si mixture. It is immedi- 
ately evident that the lowest ^^Ni-yield case, realization 3 for 
X22Ne = 0.02, is also that which is the most expanded at the 
DDT transition. 

The outcome of the deflagration appears to arise to some 
degree from the morphology of the deflagration, and there- 
fore presumably from (randomly determined) characteristics 
of the ignition condition. Cases for which the highest plume 
is near the axis, realization 1, 2 and 5, expanded less be- 
fore the detonation began than more equatorial lead plumes 
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Fig. 5. — Comparison of burning products approximately 0.1 seconds after first detonation is launched for different realization of the initial flame surface 
(columns) and for two abundances of ^-Ne (rows). Simulations are performed in axisymmetry. Fuel and burning products are indicated by color: unburned C. 
O, Ne (yellow), O-Si (red). Si-group (green), Fe-group (NSE, black). Density in g cm"^ is indicated by contours logarithmically spaced at integer powers of 10 
starting from lO' g cm"' at the outside. One extra contour (red) is added at a density of 2 X lO' g cm"'. 
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TABLE 4 

Outcomes for subset of theoretical sample 



realization 


fDDT [s] 


M(p7 > 2,?ddt) [Mq] 


Mnse [Mq] 


M56Ni [Ms] 




Mo-Si [Mq] 


energy released [lO" erg] 




1 


1.17 


1.00 


0.91 


0.79 


0.34 


0.11 


1.93 


2 


1.19 


0.96 


0.90 


0.78 


0.34 


0.11 


1.92 


3 


1.36 


0.78 


0.75 


0.64 


0.45 


0.14 


1.88 


4 


1.36 


0.77 


0.66 


0.55 


0.53 


0.15 


1.85 


5 


1.16 


0.96 


0.89 


0.77 


0.36 


0.11 


1.92 


X22Ne = 0.02 


1 


1.21 


0.87 


0.88 


0.72 


0.35 


0.11 


1.93 


2 


1.15 


0.97 


0.84 


0.68 


0.38 


0.13 


1.91 


3 


1.43 


0.54 


0.60 


0.46 


0.57 


0.17 


1.83 


4 


1.29 


0.81 


0.73 


0.58 


0.48 


0.14 


1.86 


5 


1.17 


0.92 


0.83 


0.67 


0.40 


0.13 


1.92 



SN la dependence on Ne 



15 



TABLE 5 
Initial condition amplitudes 



realization 


/ = 12 


13 


14 


15 


16 


1 


-1.20 


-1.16 


0.49 


0.21 


-1.16 


2 


-0.73 


1.17 


0.95 


-1.18 


-0.62 


3 


-0.94 


-0.64 


0.69 


0.14 


-0.61 


4 


-1.15 


-1.00 


0.03 


1.11 


0.60 


5 


0.70 


-1.69 


-1.17 


0.29 


-1.05 



as in realization 3. Realization 4 appears to form an inter- 
mediate case, which in turns makes its outcome the most 
sensitive to the inclusion of ^^Ne among these trials. This 
effect can be reasonably understood in terms of the interac- 
tion of the Rayleigh-Taylor (R-T) instability with the axisym- 
metric geometry. Plumes near the equator are dynamically 
more like rings instead of columns, therefore being more two- 
dimensional structures, and therefore are driven somewhat 
more weakly by the R-T instability. This effect may account 
for their slower rise in these simulations, and the resultant 
longer delay before detonation ignition. Additionally, equa- 
torial plumes will generally burn material more quickly due 
to a larger integrated surface and volume in axisymmetry, and 
therefore could also enhance the expansion of the star directly. 

This interaction with the geometry is important for extend- 
ing these results to three dimensions. While the direct dif- 
ferential suppression of R-T will no longer be important, the 
presence or absence of localized "spikes" in the initial condi- 
tion is likely to become the most important determining factor 
in the competition between stellar expansion and plume rise, 
and therefore the explosion yield. Thus, assumed geometry is 
no longer important, but geometric features of the initial con- 
dition serve a similar role. This also creates the possibility that 
the specific morphology imparted to the flame surface by the 
convection field in the early deflagration phase is an essential 
aspect of the outcome of the DDT scenario. Filamentary or 
sheet-like structures might act quite differently from lumps of 
burned material once the strong R-T growth phase takes over. 

5. A DEFLAGRATION DETONATION TRANSITION SUPERNOVA IN 
TWO DIMENSIONS 

In this section we describe the outcome of a typical case 
from the ensemble that will be used to study systematic ef- 
fects, the case labeled "realization 2" or "r2" for Xiij^^ = 0.02. 
This case is considered typical because the mass at high den- 
sity at the time of transition to detonation (see Section |4~2| and 
therefore the total mass of ^^Ni synthesized, O.68M0, is sim- 
ilar to the most probable values for both the ensemble of sec- 
tion |4] and, for ^^Ni m ass, the observed distribution of SNe la 
(e.g.T jHowell et al.|200 9b ). 



5.1. The Explosion 

The explosion can be roughly divided into 3 phases, the 
deflagration phase, the detonation phase, and the expansion 
phase. The deflagration phase lasts for the first 1 .2 to 1 .4 sec- 
onds of the simulation, depending on the particular ignition 
condition. During this phase a subsonic burning front, accel- 
erated by turbulence and buoyancy, burns material in the in- 
ner portion of the star and expands it appreciably. As plumes 
of burning material rise from the core, the deflagration front 
eventually reaches low enough densities to transition to a det- 
onation. The first of these to initiate a detonation will begin 
the detonation phase, which will burn the entire star within a 



few tenths of a second. Once all the material is burned, the star 
will continue to expand, beginning the expansion phase. The 
outgoing ejecta will eventually reach a state of free expansion 
in which the radial velocity of material is a linear function of 
radius and therefore the density becomes a simple function of 
time. We will now discuss each of these stages in turn for a 
typical case from our ensemble, including some details about 
how the simulation is accomplished. 

5.1.1. Deflagration 

The initial burned region and the progression of the defla- 
gration phase for realization 2 is shown in the top row of Fig- 
ure |6] The burned material is colored black throughout this 
phase, because burning in the flame at these densities always 
results in Fe-group (NSE) material within the width of the ar- 
tificial flame. Black lines show density contours on a logarith- 
mic scale in integer powers of 10 beginning at 1 in the outer 
regions and reaching to 9 in the first panel. The initial WD 
has a central density of 2.2 x 10^ g cm^. The adaptive grid 
is shown by outline of the "blocks." These are logical mesh 
regions of 16 x 16 cells that represent the smallest region that 
can be independently refined to the next level. Note also that 
PARAMESH restricts neighboring refinement regions to dif- 
fer by at most a single refinement level (a factor of 2 in cell 
size). As discussed in section [373] the background star is re- 
fined to 16 km resolution, while energy-generating regions 
are refined to 4 km. The "fluff outside the star is not refined 
except as necessary to accommodate the interior grid regions. 

5.1.2. Detonation and Unbinding 

In the simulations, the transition of the burning to detona- 
tion must be initiated artificially because thermally activated 
burning is explicitly suppressed within the RD front (the ar- 
tificially broadened flame). While section |3.1| describes a 
method for allowing thermal burning within the RD front for 
the purpose of allowing the detonation to propagate into ma- 
terial that has been "partially burned," the necessary estimate 
of temperature has so far proven to be of limited utility be- 
yond detonation propagation. Spurious detonations typically 
occur if the thermal burning proximity detection is not uti- 
lized. Detonations are initiated by setting (/)/„, the progress 
variable that represents the consumption of '^C, to 1 in one 
or several neighboring zones in one time step. Note that this 
involves no addition of energy to raise the temperature. The 
increased pressure resulting from the energy release following 
from the change in composition can then initiate a detonation 
that propagates away from the ignition point. It was found that 
for detonations at the density used here (10^ g cm"^), lighting 
a single 4x4 km cell does not always successfully launch a 
detonation. The outcome of small ignition points appeared 
to depend on the local flow characteristics, with fast-rising 
plumes being one of the more commonly problematic regions. 
This problem could be ameliorated by simply increasing the 
size of the lighted region. We found good success using a 
region with a linear radius of 8 km, resulting in the simultane- 
ous ignition of a 5 x 5 cell region, still small compared to the 
overall flow characteristics and scale for changes in density. 

For this study we have chosen to characterize the DDT point 
by a simple density criteria, pddt = 10^ g cm"-*. Whenever 
the flame front, as represented in the simulation by the RD 
front, reaches this density, we light a detonation. Because 
the top of the deflagration is generally characterized by a few 
dominant plumes, the positioning of the detonation ignition 
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Fig. 6. — Time evolution of star and burning products during tlie explosion and establishment of the expansion. 



-Ne ~ ""•^ from the sample shown in FigurelSl which is a typical case based on exploration of our sample in section|4.2| By the end of this sequence (t - 

irly in free expansion. Fuel and burning products are indicated by color: unburneo C, O, Ne {yellow), 0-Si (red). Si-group 



all burning has ceased and the star is nearly 
{green), Fe-group (NSE, black). Density in g cm 



This sequence is for realization 2 with 

3.5 s) 



is indicated by logarithmically spaced at integer powers of 10 starting from 10 g cm at the outside. One 



extra contour (red) is added at a density of 2 X 10^ g cm"''. The final panel also includes contours of radial velocity at 5, 10, 15, and 20 km s"' (thick gray). The 
adaptive mesh is indicated in the early panels by outlines of blocks of 16 X 16 cells, the smallest unit of contiguous refinement. The first detonation points are 
initiated at the moment depicted in the 5th panel at f = 1.12 s. 
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points is relatively unambiguous. While pddt quantitatively 
defines our detonation point, it is further necessary specify the 
method in which this density is used for placing the ignition 
point or points. In order to increase the chance of obtaining a 
robust detonation ignition, the ignition point is placed slightly 
outside of the rising burned region. Points are ignited when 
the RD front is approximately 64 km below the 10^ g cm"^ 
contour, at a point halfway between the RD front and the con- 
tour. Typically two points are ignited per plume, but some- 
times three if the plume is wide, as is the case for the first ig- 
nited plume in realization 2 shown in Figure [6] This method 
places the detonation points at a density of between 1 .05 and 
1.1 X lO^g cm"-'. A density contour of 2 x lO^g cm"-' is shown 
in Figure |6] in order to show that the ambiguity in detonation 
time and place, introduced by the slight difference between 
Pddt ™d the actual point of detonation ignition, should be 
a fairly small factor in the outcome of the overall explosion, 
assuming it is applied consistently. 

It should be emphasized that this treatment of detonation ig- 
nition is only a first, simplest option. The correct placement of 
the detonation ignition point is currently unknown, but at min- 
imum a more realistic condition would take into account the 
local flow characteristics in an estimate of the Gibson scale 
and compare it to the flame width. It is even possible that 
the detonation ignition is qualitatively different. For exam- 
ple, we have ignited the detonation on the "top" of the rising 
plumes, although the most pronounced mixing and tumble is 
occurring on the underside of plume caps. It may be that the 
detonation does not occur until this region reaches some den- 
sity threshold, as mixing on the plume tops is less vigorous. 
This condition would lead to a systematic delay of the deto- 
nation ignition from the methodology implemented here. We 
leave evaluation of various alternatives for ignition of the det- 
onation to future work. 

Once the detonation is ignited, most of the remainder of 
the star is burned in a few tenths of a second, between 1.12 
s and 1.6 s in realization 2. In the sequence shown in Fig- 
ure |6] additional ignition points launch from where plumes in 
the lower half-plane reach pddt- The refinement tracks the 
detonation at lower densities, but much of the higher density 
material remains refined longer because of the additional en- 
ergy release as NSE material expands and alpha particles are 
recaptured. The display of the grid has been eliminated after 
the sixth panel due to it becoming too dense in this visualiza- 
tion. Nearly all of the material in the previously deflagrated 
region is burned, mostly to Fe-group material, demonstrat- 
ing the success of the method used to allow the detonation 
to propagate into material within the RD front. Some mate- 
rial in heavily mixed regions is not fully burned, but given the 
low resolution achievable in the simulations, it is unclear how 
realistic this is. This incomplete burning does not appear to 
be a significant issue for gross yields from the explosion, but 
detailed nucleosynthesis, which will be undertaken in future 
work, will need to provide a better treatment of how deflagra- 
tion ashes interact with the detonation front. This is discussed 
more with the final yields below. 

The energetic history of the explosion is shown in Figure 
[7] where we find that, as expected, most of the energy is de- 
posited during the detonation phase, and that the star is, in 
fact, bound up to this point. Shown is the total binding energy, 
£bind, which includcs the gravitational binding energy, the in- 
ternal energy, and the kinetic energy. The exclusive source of 
energy is the nuclear energy release. The source of this energy 
is actually the change in rest mass energy of material due to 
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Fig. 7. — Dynamics of energy content during explosion. Shown are the 
total binding energy, £bmd, including the gravitational binding energy, the 
internal energy and the kinetic energy, the kinetic energy as a separate com- 
ponent, Sjjin (dashed), the accumulated energy lost in the form of neutrinos, 
£i^-loss (dotted), the change in the total rest mass energy of the star due to nu- 
clear processes, AErcst mass and the change in the total of all these energetic 
components AiJtoial- The term S^z-ioss has been multiplied by 10 for display 
purposes. 

nuclear processes. We quantify this by measuring the differ- 
ence between the rest mass of all the material on the grid and 
the same number of baryons of symmetric matter (Ye = 0.5) 
in a completely unbound state of free protons, neutrons and 
electrons: 

^rcst mass , \2/it7\ - 

--Ye{mp + me)c +{l-Ye)m„-q 



Na 



- [0.5(mp + mf) + 0.5m„] 
= (Yg - 0.5)(mp + nie - m„)c^ - q , 



(27) 



where e^^^^ mass is the rest mass energy per Na baryons (ap- 
proximately 1 g of baryons), nip, nin, and nig are the rest 
masses of the proton, neutron and electron respectively, and q 
is the average nuclear binding energy as defined in section [3T| 
The total rest mass energy is then ^rest mass = / ^rest massP'^V', 
where the integral is over the computational domain. As dis- 
cussed in section [3T| we are using the convention that density 
as represented on the grid is actually the baryon density di- 
vided by Na, and is only approximately the density in g cm"-'. 
Due to the binding energy of the unburned material, firest mass 
is in fact a large negative number at the beginning of the sim- 
ulation. Therefore, Figure [7] only displays the change during 
the simulation, AEi-est mass- It is immediately evident that this 
is the principal energy source, as its dependence is the direct 
converse of the energy deposition. 

The final contribution to the total energy balance is the en- 
ergy lost to neutrinos during the neutronization process. We 
define the cumulative neutrino loss by 



s(f) 



= / / eUt')pdVdt' , (28) 
Jo Jv 



where e^, is the neutrino loss for the local conditions as tabu- 
lated by |Seitenzahl et alT]p009 )) (see section [TT) . The energy 
loss due to neutrinos is also shown in Figure |7| where it has 
been multiplied by a factor of ten in order to be visible on this 
scale. We observe that there is no significant neutrino losses 
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during the detonation phase, indicating that the neutronization 
takes place exclusively during the deflagration. 

Adding all the energy terms together provides a useful 
check on the energy conservation of our code. The total is 
formed by Etotai = -Ebind + ^rest mass + ^tz-Ioss ■ This sum again 
has a large offset due to the binding energy of the unburned 
C-O-Ne material, so we only display the change in total en- 
ergy during the simulation. Very good energy conservation is 
observed, though there appears to be a small loss of energy 
(barely perceptible in Figure |7]) during heavy refinement and 
derefinement, apparently related to /Siest mass- This loss may be 
a unexpected interaction of interpolation with the representa- 
tion of the conserved quantities in the burning model. The loss 
is small enough to not be a concern, but will be addressed in 
the ongoing verification of our code components. 

5.1.3. Transition to Expansion 

In order to obtain an accurate evaluation of the distribution 
of ejected material in velocity, which is critical for spectral 
properties of the explosion, it is necessary to continue the 
simulation until the ejecta has reached approximate free ex- 
pansion. The detonation has propagated throughout the inte- 
rior of the star by just after f = 1 .6 s as can be seen from Figure 
|6] Some nucleosynthesis takes place later as the NSE freezes 
out, but the need for fine refinement of the energy generat- 
ing regions ends at this time since there are no propagating 
burning fronts. As a result, at 1.8 s in this simulation, the 
refinement of energy-generating regions, including regions in 
which the propagation of RD front is still taking place, is dis- 
abled. Also the refinement of non-energy-generating regions 
is coarsened from 16 to 32 km. Ideally, de-refinement with 
further expansion could take place automatically with appro- 
priate detection of the expansion of flow features. This some- 
what delicate balance was not attempted in these calculations, 
and nearly the entire ejecta remains refined at 32 km through 
the end of the calculation. 

The calculation was halted as material began to flow off of 
the grid around f = 3.5 s. The domain was chosen to be of 
such a size that this time was late enough for the ejecta to be 
in approximate free expansion. We take free expansion to be 
indicated by a linear relation between the radial velocity and 
the radius. Profiles of the radius against the radial velocity of 
the outgoing ejecta are shown in Figure [8] Shown are both 
line-out profiles at 9 = 45, 90, and 135°, as well as the aver- 
aged profile (offset by +10 Mm s"' for clarity). The average 
here is the mass weighted average radius for material mov- 
ing at a radial velocity within a bin of 0.2 km s"'. While the 
average profile has a fairly linear relation, there is some depar- 
ture from linear in the outer regions for certain directions in 
the ejecta. Additionally, the radial velocity component com- 
pletely dominates any motion at this time. Radial velocity 
contours are also shown on the final panel of Figure [6] and 
are fairly symmetric, though there is a noticeable degree of 
asymmetry in the abundance profiles. 

5.2. Explosion Yield and Ejecta Structure 

The time history of the material produced in the incinera- 
tion of the star is shown in Figure |9] In the interest of simplic- 
ity, we leave tracer particle post-processing and detailed nu- 
cleosynthesis to future work and instead just make use of the 
progress variables defined in the burning model (see Section 
|3.1| l. The ashes of the first stage of burning are predominantly 
O and Si, with some other intermediate mass elements (Ne, 
Mg), and its local mass fraction is given by Xo.si = 4>fa-4>aq- 




radial velocity (10^ km s ') 

Fig. 8. — Run of radius and mass with radial velocity alt = 3.5 s, demon- 
strating that free expansion has been approximately attained. 




time(s) 

Fig. 9. — Masses of synthesized material in time during explosion, see text 
for definitions. Note that these are estimates of final outcome, and therefore 
the relaxation of the NSE material dominantly Fe-group constituents is not 
evident in this plot. 

Further burning produces a mixture of various Si-group nu- 
clides, predominantly ^^Si, which begins in NSQE. We are 
concerned with the form that material will take in the out- 
going ejecta, after NSQE and NSE material have completely 
relaxed. As such. Figure |9] shows the final yield masses as 
simply Si- and Fe-group. The mass fraction of the Si-group 
material is given by Xsi_gi-oup = <t>aq~4>qn- Material in NSE is 
given by Xnse = 4>qn, and will consist, in the ejecta, of almost 
entirely Fe-group elements. These definitions are also used 
below to characterize the abundance profile of the ejecta and 
allow the colors displayed in figure|6]to be interpreted roughly 
in terms of the nucleosynthetic products. 
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By tracking Yg, we also have a local measure of the neu- 
tronization that has occurred, and from this we can estimate 
how much of the NSE material will be in the form of stable, 
neutron-rich Fe-group nuclides instead of ^^Ni. For = 0-5 . 
the ejected material is nearly pure ^^Ni (Meakin et al. 2008|. 
In order to estimate the local mass fraction of """Ni of a neu- 
tronized fluid element being ejected, we assume that the next 
most abundant nuclides are an equal admixture of ^''Fe and 
'''^Ni, as is observed in the tabulation of nucleosynthetic out- 
come with Ye by Meakin et al. ( 2008 1. The local mass fraction 
of ^^Ni is then estimated by 

Ye,N=[Ye-(l- 4>iiN)Ye.f] / 4>qN 

. F,.Ar-0.48212 



0.5-0.48212 







(29) 



This estimated mass fraction is also used below in discussing 
abundance profiles. 

From Figure [9] nearly all of the intermediate mass ele- 
ments are produced during the detonation phase, with a small 
amount of Si group material produced in the deflagration. 
This latter is likely to be greater for more expanded explo- 
sions with lower NSE and Ni yields. The isolation of the neu- 
tronization to the deflagration phase is also evident. While 
only half of the NSE material produced in the deflagration 
phase will be ejected as ^''Ni, nearly all of that produced in 
the detonation phase will end up as ^^Ni. Note that since 
Xiif^^ = 0.02 for this explosion, the maximum XiSNi even for 
non-neutronized material is 0.95. Although it appears all the 
burning is complete by approximately 1 .6 seconds, this is only 
true in terms of gross yields. At that time much of the Fe- 
group material will still be in NSE, and will relax as the star 
expands. The small glitch in the ^^Ni curve at 1.8 seconds is 
due to the cell-mixing upon zone de-refinement and the result- 
ing change in our estimate of the final ^^Ni yield via equation 

The elemental ejecta structure is shown by the combina- 
tion of the last panel of Figure |6] which shows radial velocity 
contours overlayed on the color-coded gross yields as defined 
above, along with Figure [TOjthat gives the spatial distribution 
of Ye at the same moment with the same radial velocity con- 
tours shown. Material below Y,, « 0.482 is expected to have 
fairly little ^*'Ni, favoring more neutron-rich nuclides instead. 
For further detailed scrutiny, we have extracted the abundance 
pattern along radial lines at several latitudes and the avera ged 
abundance profile for 0.2 km s"' bins, all shown in Figure 11 
in radial velocity. In this figure the ^*Ni mass fraction is esti- 
mated from equation (29 1. In the case of the averaged profile 



this is done before averaging. The distribution of the mass in 
radial velocity is shown in Figure |8] showing that very little 
mass is outside of approximately 15,000 km s"'. Finally the 
distribution of material in mass is shown in Figure 12 where 
the mass coordinate is defined as the mass enclosed by con- 
secutive radial velocity shells. 

5.3. Resolution Dependence 

Given the small scale of the expected flame surface struc- 
ture with respect to the grid scale in these simulations (4 km), 
and the very basic treatment of turbulent acceleration of the 
burning, some dependence on resolution is expected. The im- 
portant question is whether or not such resolution effects are 
strong enough to preclude the use of suites of 4 km simu- 
lations from being used to evaluate systematic effects. We 
find that the apparent level of uncertainty from the low reso- 
lution is acceptable when considering particular metrics, but 
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Fig. 10. — The electron molar fraction Ye att = 3.5 s after free expansion 
has been essentially attained. Contours indicated radial velocities of 5, 10, 15 
and 20 X lO' km s"' . View is the same as the last panel of Figure[6] 

that there are important dependencies on resolution that we 
hope can be resolved as treatments of subgrid flame structure 
and their usage becomes more advanced. 

The metric used in our study of systematics in section|6]be- 
low is the mass at high density at the DDT time. This time 
is defined by the moment at which the lowest density of fresh 
fuel being encountered by the flame, /O/.min, falls below pddt, 
which we are assuming here is 10^ g cm""*. It was found in 
section [43] that a good density threshold to use to define the 
mass at high density is 2 x 10^ g cm"-'. Thus it is useful to look 
at how the evolution of M(p-i > 2) as the star expands, and the 
flame rises, depends on resolution. The top panel of Figure 13 
displays the dependence of M{pj > 2) on p/.min- Evolution in 
this figure moves from right to left as the flame rises through 
the star Two different realizations from the sample developed 
in section|42]are shown. One is that described in detail in this 
section (realization 2), and the other is a more extreme case 
that produces much less NSE material (realization 3). In ad- 
dition to several resolutions, a case was run with the value of 
the parameter m that appears in the buoyancy-compensation 
prescription for the flame speed (see section [l!2| doubled. By 
simultaneously doubling this parameter and halving the reso- 
lution, the actual value of the prescribed flame speed stays the 
same. 

Our metric for measuring the expansion of the star, M{p-i > 
IjIudt), shows a fairly low sensitivity to resolution for these 
cases, especially for realization 2. There is also fairly little 
sensitivity to the value of the parameter m. This insensitiv- 
ity is likely due to the fact that the two competing processes 
at work, the pulsational expansion of the star and the rise of 
burned plumes through the star, are not overly sensitive to the 
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Fig. 1 1 . — Run of burning products with radial velocity aXt = 3.5 s. 
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Fig. 12. — Distribution of burning products with mass coordinate in ejecta. 



poorly resolved flame surface. The small-scale flame struc- 
ture is not important for the selection of the dominant plume 
because this is largely determined by the initial condition (see 
also discussion in section 4.2 1. Also, due to the delay in man- 



ifestation of the turbulence from the Rayleigh-Taylor instabil- 
ity of burned plumes, the additional flame surface area and 
concomitant energy deposition may come after the critical 
phase for the launch of the pulsational expansion of the star. 




p, . (g cm ) 

Fig. 13. — Investigation of the dependence of results on resolution. The 
top panel show s the mass above the density threshold 2 X lO' g cm"^. found 
in section [4!2] to con'elate well with total NSE yield. The curves are plot- 
ted against the minimum density of fuel being cun'ently burned by the flame, 
P/.min. such that the evolution follows the curve from right to left as the flame 
rises and the star expands. The bottom panel shows the total mass burned, 
essentially all of which is NSE material. While the value of M{p^ > 2) when 
P/.min = Pddt appears fairly robust with resolution, p rovid ing some confi- 
dence in its usage to study systematic effects in section [4~T] the total burned 
mass shows more resolution dependence. 

This delay in energy release might lessen the impact of uncer- 
tainty in the burning rate during the intermediate stages of the 
deflagration, the time when it is most difficult to model accu- 
rately. However, this situation highlights the necessity for a 
careful understanding of the very early portion of the defla- 
gration phase and the interaction of the flame surface with the 
turbulent convection field arising from the simmering phase. 

There is a stronger and more systematic dependence of the 
total burned mass on resolution, as shown by the bottom panel 
of Figure 13 However, the scatter in the consecutive resolu- 



tion cases performed makes them difficult to interpret clearly. 
In an observational sense, the resolution dependence might be 
closely related to the weakness of the neutron-enriched core 
observed in the yields from realization 2, as discussed be- 
low. The additional enhancement of burning due to turbulence 
should be most effective in the inner regions of the star where 
the intersecting wakes of rising plumes will lead to strong tur- 
bulent shearing of the flame surface. This boost may burn core 
material early enough in the expansion to enable this material 
to undergo significant electron captures. The two-dimensional 
simulations performed here (in which there is no physically 
realistic turbulence cascade) and the lack of an explicit treat- 
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ment of turbulence-flame interaction make it inappropriate for 
this to be addressed in the current study. We will simply leave 
the degree to which turbulence influences the neutronization 
of material an open question, which is not expected to depend 
on the ^^Ne content under study here due to the predominance 
of turbulence in setting the fuel consumption rate under these 
conditions. 

5.4. Brief Comparison with Observed Properties of SN la 

As seen from Figures [6j [10] [TT| and [T2j our explosion re- 
produces fairly well the observed radial velocity structure of 
the intermediate mass elements ejected in the SNe la (Fil- 
|ippenko||l997| l. We will leave a fully detailed comparison 
until full abundances can be obtained from post-processing 
tracer particles, but many important features warrant men- 
tioning here. There is a prominence of Si in abundance for 
the velocity range of roughly 10 to 15 Mm s~', accounting for 
nearly QAMq of the final ejecta, with some Si-group material 
extending down as far as 8 Mm s"'. This is broadly consis- 
tent with the spectral evolution of SNe la, though possibly 
1-2 Mm s"' high in velocity compared to the spectral evolu- 
tion of the most standard cases like 1994D ( |Patat et al.|1996| l. 
Since the mass distribution (Figure |8]l is weighted toward 
lower velocities, conclusive comparison will require radiative 
transfer calculations. Also, this is only a single case of many 
possible, which may have Si at lower velocities. The ejecta is 
fairly symmetric, with a slight asymmetry arising from the lo- 
cation and timing of the detonation sites. While the degree of 
asymmetry is likely not too unrealistic, the pattern imprinted 
may be a bit different in a three-dimensional treatment. It does 
appear that the velocity extent of Si features near the photo- 
sphere would be viewing-angle-dependent. Note that the jet- 
like columns of Fe-group material on the symmetry axes are 
unrealistic, but appear to have a fairly small contribution to the 
averaged profile shown in the last panel of Figure 1 1 Detailed 



radiative transfer calculations will be required to determine if 
the patchiness of the border between the Si and Fe-group ma- 
terial provides a good match to polarimetric observations (e.g. 
Wang et al.|2003] l. 



The distribution of neutron-rich material demonstrated in 
Figures[TO]and[TT]compares less well with the distribution in- 
ferred from observations ( jStehle et al. 2005, M azzah et al.| 
2008| l. On the positive side, the impurity created by the 



distribution of deflagration material at velocities around 3- 
8 Mm s~' dilutes the pure ^^Ni to a degree that provides a 
reasonable agreement with observations. However, notably 
absent is the predominance of stable Fe at low velocities, 
< 3 Mm s~', which is inferred from observations of 2002bo 
( IStehle et al.|2005| ) and 2004eo ( |Mazzah et al.|2008| l. While 
the treatment of the ejecta evolution leading to this inferred 
abundance is approximate, it seems unlikely that such a stark 
difference could be reconciled without a change in the model. 
The most apparent reason for this shortcoming of the model 
is the flame treatment being utilized here. The buoyancy- 
compensation form of |Khokhlov| ( |1995] l, which we are using, 
does not account for interaction of the turbulent wakes with 
other burning fronts. Thus as plumes rise out of the core, the 
turbulence they leave behind does not lead to the enhance- 
ment of flame spreading that it should. This shortcoming will 
be addressed in future work, with more sophisticated treat- 
ments of the turbulence-flame interaction ( jColin et al.||2000t 
ISchmidt et al. 2006). Investigation of this effect should be 
undertaken with attention to the detailed nuclear products, as 
extensive central burning early in the deflagration phase can 



lead to "over"-neutronization (e.g. Brachwitz et al.|2000 1. 

Many SNe la have been o bserved with high-velocity Ca 
features ( Mazzali et al.|2003] l. This material is generally pro- 
duced during the transition from NSQE to NSE (e.g. Nomoto] 
|et al.| |1984) and thus should appear at the border between 
Si-group (green) and NSE (black) regions in Figure [6] In 
this simulation, these yields are restricted to velocities < 
13,000 km s"', where the main component also lies in ob- 
served supernovae. Some features, however, do indicate that 
this material might be produced at higher velocities by inter- 
secting detonation fronts. The case of the symmetry axes, 
while an artifact here, suggests that points where 3 detona- 
tions spreading through the star meet might be good candidate 
sites for the formation of such high-velocity features. This is 
reinforced by features seen in some detonation interactions in 
other realizations (see the comparison of ignition conditions 
in Figure |5]). Realistic tests will await DDT calculations in 
three dimensions. Additionally, we have assumed that deto- 
nations are ignited when the first flame front passes through 
the DDT density. This may be too stringent a condition, and 
some number of deflagration plumes may extend well above 
this transition region (which also corresponds closely to the 
division between Si- and NSE-rich ashes) before the detona- 
tion takes place. 

6. SYSTEMATIC EFFECTS OF -^NE ARISING FROM DYNAMICS 

We have implemented the current study as an attempt to 
determine whether the presence of ^^Ne has a significant sys- 
tematic effect on the expansion state of the star at the time of 
first detonation. In light of the relation between the expan- 
sion and final yields discussed in Section 4.2 this is a likely 
candidate for an effect of ^^Ne that can compete with that ex- 
pected due to the simple overall neutron excess highlighted by 
[Timmes et al.| ( [2003) l. Our focus on the degree of expansion, 
quantified by Mip-j > 2, foDx), in such specificity is motivated 
by a desire to separate the contribution of ^-Ne to the overall 
nucleosynthesis into a dynamical component distinguishable 
from the gross neutron excess. The use of the outcome of 
the deflagration phase alone also allows many more cases to 
be run at the same computational cost. This is very impor- 
tant given the high variance in the outcome of the supernovae. 
Averaging over many realizations is necessary to obtain un- 
certainties that are small compared to the expected effect of 
neutron excess. Additionally, targeting our study so specifi- 
cally gives more power for tertiary conclusions about the im- 
pact of other effects that might modify the propagation speed 
of the burning front. 

Variation due to ^^Ne was explored by simulating the de- 
flagration phase for 20 realizations with Xii^e = 0, 0.02, and 
0.04. As discussed in section |2] the carbon abundance was 
kept fixed at 0.5 by mass as was the central density of the ini- 
tial WD. This results in the WD being slightly less massive 
with a higher ^^Ne content. Our WD, as mapped on the ini- 
tial grid, is 2.7417, 2.7314 or 2.7211 x 10^^ g for Xn^, = 0, 
0.02, and 0.04 respectively. Samples with different ^^Ne con- 
tent are compared in Figure [14] where M{pj > 2,fDDT) for 
the ^^Ne enhanced cases is plotted against the same outcome 
for no ^^Ne. The dashed line indicates a 1-1 relation. Sys- 
tematic effects should arise as a departure from this relation 
on average. The two nearly coincident points with error bars 
indicate the averages of the samples with different -^Ne abun- 
dances. The outer error bars indicate the standard deviation 
and the inner the standard deviation of the mean. The sample 
with X22Ne = 0.04 has a slightly higher average of 0.85 IM© 
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Fig. 14. — Relation of degree of pre-expansion for cases with X22fv[e = 0.02 
(circles) and 0.04 (squares), to a case with no ^^Ne. Degree of pre-expansion 
is indicated by the mass above a density of 2 X lO' g cm"' when the first det- 
onation initiates, which is an indicator of how much Fe-group material will 
be ejected in the supernova. The point with vertical and horizontal error bars 
indicates the average, sample standard deviation (outer eiTor bars) and stan- 
dard deviation of the mean (inner error bars). The sample with Xizjvje = " "4 
has a slightly higher average. There is not a significant difference between 
the average expansion in the three cases. 

compared to O.846M0 for Xn-^^ = 0.02, and both of these are 
somewhat lower than the 0.88 IM0 for no ^^Ne. Standard de- 
viations of the mean are O.O42M0 for both of the nonzero 
^^Ne cases and O.O46M0 for the zero ^^Ne. 

The overlap of the inner error bars with the 1-1 relation in- 
dicates that there is no statistically significant impact of ^^Ne 
content on the expansion prior to DDT over this range of 
mass fractions. On the basis of neutron excess (ITimmes et al.| 
2003 I, the decrease in the ^^Ni mass between Ne fractions 
of and 0.04 is ~ 10%. Our results preclude a difference 
in the overall NSE production that could compete with this 
change, finding it to be < 5%. The more physically moti- 
vated interval between 0.02, appropriate for a zero metallicity 
progenitor, and 0.04, appropriate for an approximately solar 
metallicity progenitor ( Piro & Bildsten|2 008 ; Chamulak et al. 
2008| l, shows a very small difference. Note that changes in 
X22[,je by 0.02 can lead to changes in individual cases by as 
much as O.2M0, but that this effect is washed out by aver- 
aging over a variety of ignition conditions. This implies that 
the important controlling processes are large-scale plume rise 
and stellar expansion, and that these are not influenced in a 
systematic way by modest changes in the flame propagation 
like those introduced by -^Ne. This, in fact, bodes well for 
the robustness of simulations in the DDT scenario, as many 
features of the explosion don't depend on accurate treatment 
of physical phenomena at all scales. However, far more inves- 
tigation is necessary before such a robustness can be actually 
claimed. 

Returning to the cases shown in Figure |5] the addition of 
^^Ne does not change the major morphology of the deflagra- 
tion. Generally the same dominant plumes are observed at 
both ^^Ne fractions. However, the substructure of the ris- 
ing plumes can undergo some fairly extensive modification, 
which can also influence the timing of their rise and therefore 
the detonation. The modification likely arises from influence 
of the flame speed on the early evolution of structure on the 
smallest scales of the initial flame surface. The absence of 



a systematic effect implies that such modification is equally 
effective at suppressing dominant plumes and enhancing sub- 
dominant ones. This supposition is demonstrated anecdotally 
in Figure [5] where, for example, in realization 5, enhancing 
the ^^Ne enhances one secondary plume and suppresses an- 
other. This effect might be a manifestation of the tendency of 
R-T to cause perturbations to grow at the expense of others, 
such that selection of the dominant perturbation is less impor- 
tant than the dynamics of the dominant perturbation once one 
"wins". 

There is still an important way for ^^Ne to influence the 
dynamics of the explosion not addressed here. This study 
treated a sample of ignition conditions with fixed harmonic 
content. If, however, the change in the flame speed due to 
^^Ne has a significant impact during the earliest stages of the 
deflagration, it may change the effective harmonic content of 
an abstracted initial condition as implemented here. However, 
the apparent lack of sensitivity, on average, to changes in the 
small scale structure of deflagration morphology argue against 
this being important. Thus, the important features of the ig- 
nition condition may be set by the large scale flows in the 
convective core and the relative position of the initial spark. 
This subject deserves attention as simulations of the early part 
of the deflagration phase are undertaken. 

Although we have held other parameters fixed for this study, 
as discussed in section|2] the ^^Ne composition is not the only 
parameter which may vary with the metallicity of the parent 
stellar population. Notably the central carbon fraction and ig- 
nition density are expected to depend on both the metallicity 
and the main-sequence mass of the progenitor star. The igni- 
tion density may additionally depend on the accretion history 
of the progenitor WD. A full accounting of the systematic de- 
pendence of theoretical SNe la will therefore have to await 
evaluation of these further parameters, enabled by the frame- 
work presented here. 

7. CONCLUSIONS AND FUTURE WORK 

We have extended a deflagration ignition condition of con- 
strained asymmetry, with no harmonic components less than 
^ = 12 in the initial flame surface, to develop a well-defined 
random sample, and use this to construct a basic framework 
for formal study of systematic effects in SN la. Our theoret- 
ical sample has a ^^Ni mass average and standard deviation 
of 0.7 and O.IIM0 (0.62 and O.IIM0) for a ^^Ne fraction of 
(0.02) by mass, and ^''Ni yield ranging between 0.45 and 
O.8M0. This sample compares well with that of observed 
SN la (e.g. [Howell et al. 2009b ), and should be somewhat tun- 
able by varying the harmonic content of the initial condition 
and/or the DDT density. Our framework uses this sample to 
enable statistically well-defined studies of both physical and 
theoretical parameters of the SN la explosion simulation. A 
great virtue of this methodology is that it includes the inher- 
ent diversity of SN la outcomes and the real-world challenges 
that brings for the study of systematic effects. The frame- 
work is expected to undergo refinement and extension to make 
the simulations of the SN la increasingly realistic. Important 
studies for these improvements include improved turbulence- 
flame interaction models and a more careful characterization 
of the relation between spectral content in the initial condition 
and simulation outcome. 

The outcome of a typical two-dimensional simulation of an 
SN la from our sample with a DDT density of 10' g cm"-^ 
was presented in detail and compared to observations. This 
model explosion provides a reasonable reproduction of many 
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features of observed SN la. Intermediate mass elements are 
dominant at velocities above 10 Mm s~', and a mixture of 
neutron-enriched Fe-group and ""^Ni dominates at lower ve- 
locities. Material is fairly well constrained to layers in ve- 
locity, with a modest degree of overall asymmetry. There is 
also some large-scale mixing at the border between Si-group 
and Fe-group dominated material in the velocity range ~ 8- 
1 1 Mm s"' , due to tops of the highest plumes during the de- 
flagration phase. For this study we have relied upon the pa- 
rameterized burning model used in the hydrodynamic simu- 
lation in order to give a gross measure of the nucleosynthetic 
outcome. Future work will proceed with post-processing La- 
grangian particle tracks, yielding full nucleosynthetic infor- 
mation. 

Besides the many features that match observations well, 
it appears that our simple prescription for turbulence-flame 
interaction, which only compensates for adverse effects of 
buoyancy on the integrity of the burning front, causes an un- 
derproduction of stable Fe-group elements at low velocities 
(< 3 Mm s~'). Improved treatments that have been proposed 
in the literature (Col in et al. 2000 , Schmidt et al. 200 6) will 
be included in future work, allowing a direct comparison. No- 
tably, however, we do not find a strong resolution dependence 
for the competition between the stellar expansion and plume 
rise that sets the star's density structure at detonation and, 
therefore, the overall yield of NSE and intermediate mass ma- 
terial. This result indicates that our prescription is sufficient 
for its main intended purpose of allowing an approximate cap- 
ture of the dynamics of large-scale plumes. 

Although we briefly reviewed the relevant dependencies 
that must be considered in order to evaluate the full systematic 
effect of stellar population metallicity (Section[2]i, for the first 
application of the framework we focus on a question of more 
well-defined scope. We fix the carbon content of the white 
dwarf to half by mass, and assume the same accretion history, 
so that the ignition density is the same. We also assume the 
DDT density has a single value, 10^ g cm"-'. These assump- 
tions make our study specifically tailored to understand how 
dynamical effects might change the simple relations hip be- 



tween ^^Ne content and ^^Ni production presented by |Timmes 
[etar] ( [2003l l. 

We proceed by calibrating an indicator of the degree of ex- 
pansion that takes place prior to the detonation ignition, which 
is chosen to be representative of the total amount of NSE 
material produced in the explosion. We find that the mass 
above a density of 2 x 10^ g cm"^ when the DDT takes place, 
M(pj > 2,fDDT), provides a good indicator of the final NSE 
mass. This both provides a very direct probe of the dynamical 
contribution of ^^Ne and can be evaluated efficiently for many 
realizations in the theoretical sample, giving good statistics. 
By averaging over 20 realizations of initial condition, we 
find that M(p^ > 2,fDDT) is 0.881, 0.846 and 0.851 Mq for 
X22^g = 0, 0.02, and 0.04 respectively all with a standard de- 
viations of the mean of approximately O.O4M0. These results 
are shown in Figure 15 as estimated NSE masses (squares). 
There is no statistically significant dependence of the star's 
expansion, and therefore the total NSE mass produced, on the 
^^Ne content. In individual cases, however, M{pj > 2,fDDT) 
can vary by as much as O.2M0 due to a change in Xiij^^ of 
0.02. Further, any possible effect is smaller t han the 10% re- 
duction in ''^Ni mass predicted by Timmes et al. (2003) on the 
basis of neutron excess for this range of ^^Ne fraction assum- 
ing constant NSE mass. Their prediction, which assumes that 
NSE mass is independent of ^^Ne content, is shown by dashed 




Ne mass fraction 



Fig. 15.- 



Comparison of our simulation results and the trends pre- 
dicted trom vaiiation in inheiited neutron excess alone. jTimmes et al.|{2003) 
predicts that for constant NSE mass produced in the SN la (upper aashed 
line), the ^*Ni mass produced decreases linearly with progenitor ^^Ne content 
(lower dashed line). Our results from 20 simulations of the SN la deflagration 
phase at each abundance indicate that the average NSE masses (squares) ai'e 
consistent with being independent of -^Ne content, an d that the trend of ^^N i 
(circles) is therefore consistent with that predicted by |Timmes et al.H2003) . 
Only a subset of simulations were run far enough to obtain final ■'"Ni yields. 



lines in Figure 



15 



Also shown is the average ''^Ni mass yields 
from the 5 simulations at each of Xn-f^^ = and 0.02. For these, 
the error bars are obtained using the standard deviation of the 
full sample, as this small subset underestimates the variance. 
These are consistent with the reduction due to neutron excess, 
although the statistical uncertainty is fairly large. 

Our studies point to the morphology of the ignition con- 
dition as being the dominant dynamical driver of the ^''Ni 
yield of the explosion. In two dimensions this is manifested 
by whether the flow is dominated by an equatorial or polar 
plume, but there should be a analogous morphological criteria 
in three dimensions. This points to the importance of the very 
early deflagration phase, during which the region taken here 
as already burned is formed. The interaction of flame prop- 
agation and turbulence in this region will set the harmonic 
content of the initial condition for a study such as this one. 
This provides a further opportunity for ^^Ne to be important 
through its effect on the flame propagation speed. 

Finally it is important to keep in mind that the effect of 
^^Ne on both the propagation speed and width of the flame is 
expected to change the DDT density ( Chamulak et al.|2008[ l 



even if the dynamics is not affected on average. This will be 
addressed in future work. 
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